program QR
  implicit none
  INTEGER MMAX,NMAX, LDA, LWORK,nb
  PARAMETER (MMAX=3,NB=64,NMAX=3, lda= Mmax)
  PARAMETER (LWORK=2*NMAX+(NMAX+1)*NB) !NMAX=NS
  INTEGER I, IFAIL, INFO, J, N,  M
  parameter (N=3)
  DOUBLE PRECISION A(Lda,NMAX), WORK(LWORK), TAU(NMAX)
  INTEGER JPVT(NMAX)
  JPVT = 0
  A(1,1)=1.0
  A(1,2)=4.0
  A(1,3)=7.0
  A(2,1)=2.0
  A(2,2)=5.0
  A(2,3)=8.0
  A(3,1)=3.0
  A(3,2)=6.0
  A(3,3)=9.0
  CALL DGEQP3(N, N, A, LDA,jpvt,  TAU, WORK, LWORK,INFO) !qr decomposition
  write(*,99999) A
  CALL DORGQR(N, N, N, A, LDA, TAU, WORK, LWORK,INFO)    !generate Q
  write(*,99999) A
  99999 FORMAT (1X,(3f10.5))
  end
结果是
  -13.92839   0.38226   0.43004
  -3.58979  -1.05518   0.78897
  -8.75909  -0.52759   0.00000

   -0.50257  -0.57437  -0.64616
   0.76207   0.05862  -0.64483
   0.40825  -0.81650   0.40825
那按我的理解第一个矩阵的上三角为R矩阵，下面的矩阵就是Q，但是为什么Q与R的乘积不是输入矩阵A？大家帮忙看看，谢谢啦！