Télécharger la plaquette Cast3M


fr en

Télécharger qrsolv.eso

Retour à la liste

Numérotation des lignes : 

qrsolv

  1. C QRSOLV    SOURCE    PV090527  23/01/30    21:15:04     11574          
  2.       subroutine daxpy ( n, da, dx, incx, dy, incy )
  3.  
  4. c*****************************************************************************80
  5. c
  6. c! DAXPY computes constant times a vector plus a vector.
  7. c
  8. c  Discussion:
  9. c
  10. c    This routine uses double precision real arithmetic.
  11. c
  12. c    This routine uses unrolled loops for increments equal to one.
  13. c
  14. c  Licensing:
  15. c
  16. c    This code is distributed under the GNU LGPL license. 
  17. c
  18. c  Modified:
  19. c
  20. c    16 May 2005
  21. c
  22. c  Author:
  23. c
  24. c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
  25. c    David Kincaid, Fred Krogh.
  26. c    FORTRAN90 version by John Burkardt.
  27. c
  28. c  Reference:
  29. c
  30. c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
  31. c    LINPACK User's Guide,
  32. c    SIAM, 1979,
  33. c    ISBN13: 978-0-898711-72-1,
  34. c    LC: QA214.L56.
  35. c
  36. c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
  37. c    Algorithm 539, 
  38. c    Basic Linear Algebra Subprograms for Fortran Usage,
  39. c    ACM Transactions on Mathematical Software, 
  40. c    Volume 5, Number 3, September 1979, pages 308-323.
  41. c
  42. c  Parameters:
  43. c
  44. c    Input, integer N, the number of elements in DX and DY.
  45. c
  46. c    Input, real*8 DA, the multiplier of DX.
  47. c
  48. c    Input, real*8 DX(*), the first vector.
  49. c
  50. c    Input, integer INCX, the increment between successive 
  51. c    entries of DX.
  52. c
  53. c    Input/output, real*8 DY(*), the second vector.
  54. c    On output, DY(*) has been replaced by DY(*) + DA * DX(*).
  55. c
  56. c    Input, integer INCY, the increment between successive 
  57. c    entries of DY.
  58.  
  59.         implicit real*8 (a-h,o-z)
  60.         implicit integer (i-n)
  61.  
  62.         real*8 da
  63.         real*8 dx(*)
  64.         real*8 dy(*)
  65.         integer i
  66.         integer incx
  67.         integer incy
  68.         integer ix
  69.         integer iy
  70.         integer m
  71.         integer n
  72.  
  73.         if ( n .le. 0 ) then
  74.           return
  75.         end if
  76.  
  77.         if ( da .eq. 0.0D+00 ) then
  78.           return
  79.         end if
  80.  
  81. c  Code for unequal increments or equal increments
  82. c  not equal to 1.
  83.  
  84.         if ( incx /= 1 .or. incy /= 1 ) then
  85.  
  86.           if ( 0 .le. incx ) then
  87.             ix = 1
  88.           else
  89.             ix = ( - n + 1 ) * incx + 1
  90.           end if
  91.  
  92.           if ( 0 .le. incy ) then
  93.             iy = 1
  94.           else
  95.             iy = ( - n + 1 ) * incy + 1
  96.           end if
  97.  
  98.           do i = 1, n
  99.             dy(iy) = dy(iy) + da * dx(ix)
  100.             ix = ix + incx
  101.             iy = iy + incy
  102.           end do
  103.  
  104. c  Code for both increments equal to 1.
  105.  
  106.         else
  107.  
  108.           m = mod ( n, 4 )
  109.  
  110.           dy(1:m) = dy(1:m) + da * dx(1:m)
  111.  
  112.           do i = m+1, n, 4
  113.             dy(i  ) = dy(i  ) + da * dx(i  )
  114.             dy(i+1) = dy(i+1) + da * dx(i+1)
  115.             dy(i+2) = dy(i+2) + da * dx(i+2)
  116.             dy(i+3) = dy(i+3) + da * dx(i+3)
  117.           end do
  118.  
  119.         end if
  120.  
  121.         return
  122.       end
  123.       function dnrm2 ( n, x, incx )
  124.  
  125. c*****************************************************************************80
  126. c
  127. c! DNRM2 returns the euclidean norm of a vector.
  128. c
  129. c  Discussion:
  130. c
  131. c    This routine uses double precision real arithmetic.
  132. c
  133. c     DNRM2 ( X ) = sqrt ( X' * X )
  134. c
  135. c  Licensing:
  136. c
  137. c    This code is distributed under the GNU LGPL license. 
  138. c
  139. c  Modified:
  140. c
  141. c    16 May 2005
  142. c
  143. c  Author:
  144. c
  145. c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
  146. c    David Kincaid, Fred Krogh.
  147. c    FORTRAN90 version by John Burkardt.
  148. c
  149. c  Reference:
  150. c
  151. c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
  152. c    LINPACK User's Guide,
  153. c    SIAM, 1979,
  154. c    ISBN13: 978-0-898711-72-1,
  155. c    LC: QA214.L56.
  156. c
  157. c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
  158. c    Algorithm 539, 
  159. c    Basic Linear Algebra Subprograms for Fortran Usage,
  160. c    ACM Transactions on Mathematical Software,
  161. c    Volume 5, Number 3, September 1979, pages 308-323.
  162. c
  163. c  Parameters:
  164. c
  165. c    Input, integer N, the number of entries in the vector.
  166. c
  167. c    Input, real*8 X(*), the vector whose norm is to be computed.
  168. c
  169. c    Input, integer INCX, the increment between successive 
  170. c    entries of X.
  171. c
  172. c    Output, real*8 DNRM2, the Euclidean norm of X.
  173. c
  174.         implicit real*8 (a-h,o-z)
  175.         implicit integer (i-n)
  176.  
  177.         real*8 absxi
  178.         real*8 dnrm2
  179.         integer incx
  180.         integer ix
  181.         integer n
  182.         real*8 norm
  183.         real*8 scale
  184.         real*8 ssq
  185.         real*8 x(*)
  186.  
  187.         if ( n .lt. 1 .or. incx .lt. 1 ) then
  188.  
  189.           norm  = 0.0D+00
  190.  
  191.         else if ( n .eq. 1 ) then
  192.  
  193.           norm  = abs ( x(1) )
  194.  
  195.         else
  196.  
  197.           scale = 0.0D+00
  198.           ssq = 1.0D+00
  199.  
  200.           do ix = 1, 1 + ( n - 1 ) * incx, incx
  201.             if ( x(ix) /= 0.0D+00 ) then
  202.               absxi = abs ( x(ix) )
  203.               if ( scale .lt. absxi ) then
  204.                 ssq = 1.0D+00 + ssq * ( scale / absxi )**2
  205.                 scale = absxi
  206.               else
  207.                 ssq = ssq + ( absxi / scale )**2
  208.               end if
  209.             end if
  210.           end do
  211.           norm  = scale * sqrt ( ssq )
  212.         end if
  213.  
  214.         dnrm2 = norm
  215.  
  216.         return
  217.       end
  218.       subroutine dqrank ( a, lda, m, n, tol, kr, jpvt, qraux, work )
  219.  
  220. c*****************************************************************************80
  221. c
  222. c! DQRANK computes the QR factorization of a rectangular matrix.
  223. c
  224. c  Discussion:
  225. c
  226. c    This routine is used in conjunction with sqrlss to solve
  227. c    overdetermined, underdetermined and singular linear systems
  228. c    in a least squares sense.
  229. c
  230. c    DQRANK uses the LINPACK subroutine DQRDC to compute the QR
  231. c    factorization, with column pivoting, of an M by N matrix A.
  232. c    The numerical rank is determined using the tolerance TOL.
  233. c
  234. c    Note that on output, ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
  235. c    of the condition number of the matrix of independent columns,
  236. c    and of R.  This estimate will be .le. 1/TOL.
  237. c
  238. c  Reference:
  239. c
  240. c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
  241. c    LINPACK User's Guide,
  242. c    SIAM, 1979,
  243. c    ISBN13: 978-0-898711-72-1,
  244. c    LC: QA214.L56.
  245. c
  246. c  Parameters:
  247. c
  248. c    Input/output, real*8 A(LDA,N).  On input, the matrix whose
  249. c    decomposition is to be computed.  On output, the information from DQRDC.
  250. c    The triangular matrix R of the QR factorization is contained in the
  251. c    upper triangle and information needed to recover the orthogonal
  252. c    matrix Q is stored below the diagonal in A and in the vector QRAUX.
  253. c
  254. c    Input, integer LDA, the leading dimension of A, which must
  255. c    be at least M.
  256. c
  257. c    Input, integer M, the number of rows of A.
  258. c
  259. c    Input, integer N, the number of columns of A.
  260. c
  261. c    Input, real*8 TOL, a relative tolerance used to determine the
  262. c    numerical rank.  The problem should be scaled so that all the elements
  263. c    of A have roughly the same absolute accuracy, EPS.  Then a reasonable
  264. c    value for TOL is roughly EPS divided by the magnitude of the largest
  265. c    element.
  266. c
  267. c    Output, integer KR, the numerical rank.
  268. c
  269. c    Output, integer JPVT(N), the pivot information from DQRDC.
  270. c    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
  271. c    independent to within the tolerance TOL and the remaining columns
  272. c    are linearly dependent.
  273. c
  274. c    Output, real*8 QRAUX(N), will contain extra information defining
  275. c    the QR factorization.
  276. c
  277. c    Workspace, real*8 WORK(N).
  278. c
  279.         implicit real*8 (a-h,o-z)
  280.         implicit integer (i-n)
  281.  
  282.         integer lda
  283.         integer n
  284.  
  285.         real*8 a(lda,n)
  286.         integer j
  287.         integer job
  288.         integer jpvt(n)
  289.         integer k
  290.         integer kr
  291.         integer m
  292.         real*8 qraux(n)
  293.         real*8 tol
  294.         real*8 work(n)
  295.  
  296.         jpvt(1:n) = 0
  297.         job = 1
  298.  
  299.         call dqrdc ( a, lda, m, n, qraux, jpvt, work, job )
  300.  
  301.         kr = 0
  302.         k = min ( m, n )
  303.  
  304.         do j = 1, k
  305.           if ( abs ( a(j,j) ) .le. tol * abs ( a(1,1) ) ) then
  306.             return
  307.           end if
  308.           kr = j
  309.         end do
  310.  
  311.         return
  312.       end
  313.       subroutine dqrdc ( a, lda, n, p, qraux, jpvt, work, job )
  314.  
  315. c*****************************************************************************80
  316. c
  317. c! DQRDC computes the QR factorization of a real rectangular matrix.
  318. c
  319. c  Discussion:
  320. c
  321. c    DQRDC uses Householder transformations.
  322. c
  323. c    Column pivoting based on the 2-norms of the reduced columns may be
  324. c    performed at the user's option.
  325. c
  326. c  Reference:
  327. c
  328. c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
  329. c    LINPACK User's Guide,
  330. c    SIAM, 1979,
  331. c    ISBN13: 978-0-898711-72-1,
  332. c    LC: QA214.L56.
  333. c
  334. c  Parameters:
  335. c
  336. c    Input/output, real*8 A(LDA,P).  On input, the N by P matrix
  337. c    whose decomposition is to be computed.  On output, A contains in
  338. c    its upper triangle the upper triangular matrix R of the QR
  339. c    factorization.  Below its diagonal A contains information from
  340. c    which the orthogonal part of the decomposition can be recovered.
  341. c    Note that if pivoting has been requested, the decomposition is not that
  342. c    of the original matrix A but that of A with its columns permuted
  343. c    as described by JPVT.
  344. c
  345. c    Input, integer LDA, the leading dimension of the array A.
  346. c    LDA must be at least N.
  347. c
  348. c    Input, integer N, the number of rows of the matrix A.
  349. c
  350. c    Input, integer P, the number of columns of the matrix A.
  351. c
  352. c    Output, real*8 QRAUX(P), contains further information required
  353. c    to recover the orthogonal part of the decomposition.
  354. c
  355. c    Input/output, integer JPVT(P).  On input, JPVT contains
  356. c    integers that control the selection of the pivot columns.  The K-th
  357. c    column A(*,K) of A is placed in one of three classes according to the
  358. c    value of JPVT(K).
  359. c      .gt. 0, then A(K) is an initial column.
  360. c      = 0, then A(K) is a free column.
  361. c      .lt. 0, then A(K) is a final column.
  362. c    Before the decomposition is computed, initial columns are moved to
  363. c    the beginning of the array A and final columns to the end.  Both
  364. c    initial and final columns are frozen in place during the computation
  365. c    and only free columns are moved.  At the K-th stage of the
  366. c    reduction, if A(*,K) is occupied by a free column it is interchanged
  367. c    with the free column of largest reduced norm.  JPVT is not referenced
  368. c    if JOB .eq. 0.  On output, JPVT(K) contains the index of the column of the
  369. c    original matrix that has been interchanged into the K-th column, if
  370. c    pivoting was requested.
  371. c
  372. c    Workspace, real*8 WORK(P).  WORK is not referenced if JOB .eq. 0.
  373. c
  374. c    Input, integer JOB, initiates column pivoting.
  375. c    0, no pivoting is done.
  376. c    nonzero, pivoting is done.
  377. c
  378.         implicit real*8 (a-h,o-z)
  379.         implicit integer (i-n)
  380.  
  381.         integer lda
  382.         integer n
  383.         integer p
  384.  
  385.         real*8 a(lda,p)
  386.         integer jpvt(p)
  387.         real*8 qraux(p)
  388.         real*8 work(p)
  389.         integer j
  390.         integer job
  391.         integer jp
  392.         integer l
  393.         integer lup
  394.         integer maxj
  395.         real*8 maxnrm
  396.         real*8 nrmxl
  397.         integer pl
  398.         integer pu
  399.         real*8 ddot
  400.         real*8 dnrm2
  401.         logical swapj
  402.         real*8 t
  403.         real*8 tt
  404.  
  405.         pl = 1
  406.         pu = 0
  407.  
  408. c  If pivoting is requested, rearrange the columns.
  409.  
  410.         if ( job /= 0 ) then
  411.  
  412.           do j = 1, p
  413.  
  414.             swapj = 0 .lt. jpvt(j)
  415.  
  416.             if ( jpvt(j) .lt. 0 ) then
  417.               jpvt(j) = - j
  418.             else
  419.               jpvt(j) = j
  420.             end if
  421.  
  422.             if ( swapj ) then
  423.  
  424.               if ( j /= pl ) then
  425.                 call dswap ( n, a(1,pl), 1, a(1,j), 1 )
  426.               end if
  427.  
  428.               jpvt(j) = jpvt(pl)
  429.               jpvt(pl) = j
  430.               pl = pl + 1
  431.  
  432.             end if
  433.  
  434.           end do
  435.  
  436.           pu = p
  437.  
  438.           do j = p, 1, -1
  439.  
  440.             if ( jpvt(j) .lt. 0 ) then
  441.  
  442.               jpvt(j) = - jpvt(j)
  443.  
  444.               if ( j /= pu ) then
  445.                 call dswap ( n, a(1,pu), 1, a(1,j), 1 )
  446.                 jp = jpvt(pu)
  447.                 jpvt(pu) = jpvt(j)
  448.                 jpvt(j) = jp
  449.               end if
  450.  
  451.               pu = pu - 1
  452.  
  453.             end if
  454.  
  455.           end do
  456.  
  457.         end if
  458.  
  459. c  Compute the norms of the free columns.
  460.  
  461.         do j = pl, pu
  462.           qraux(j) = dnrm2 ( n, a(1,j), 1 )
  463.         end do
  464.  
  465.         work(pl:pu) = qraux(pl:pu)
  466.  
  467. c  Perform the Householder reduction of A.
  468.  
  469.         lup = min ( n, p )
  470.  
  471.         do l = 1, lup
  472.  
  473. c  Bring the column of largest norm into the pivot position.
  474.  
  475.           if ( (pl .le. l) .and. (l .lt. pu) ) then
  476.  
  477.             maxnrm = 0.0D+00
  478.             maxj = l
  479.             do j = l, pu
  480.               if ( maxnrm .lt. qraux(j) ) then
  481.                 maxnrm = qraux(j)
  482.                 maxj = j
  483.               end if
  484.             end do
  485.  
  486.             if ( maxj /= l ) then
  487.               call dswap ( n, a(1,l), 1, a(1,maxj), 1 )
  488.               qraux(maxj) = qraux(l)
  489.               work(maxj) = work(l)
  490.               jp = jpvt(maxj)
  491.               jpvt(maxj) = jpvt(l)
  492.               jpvt(l) = jp
  493.             end if
  494.  
  495.           end if
  496.  
  497. c  Compute the Householder transformation for column L.
  498.  
  499.           qraux(l) = 0.0D+00
  500.  
  501.           if ( l /= n ) then
  502.  
  503.             nrmxl = dnrm2 ( n-l+1, a(l,l), 1 )
  504.  
  505.             if ( nrmxl /= 0.0D+00 ) then
  506.  
  507.               if ( a(l,l) /= 0.0D+00 ) then
  508.                 nrmxl = sign ( nrmxl, a(l,l) )
  509.               end if
  510.  
  511.               call dscal ( n-l+1, 1.0D+00 / nrmxl, a(l,l), 1 )
  512.               a(l,l) = 1.0D+00 + a(l,l)
  513.  
  514. c  Apply the transformation to the remaining columns, updating the norms.
  515.  
  516.               do j = l + 1, p
  517.  
  518.                 t = - ddot ( n-l+1, a(l,l), 1, a(l,j), 1 ) / a(l,l)
  519.                 call daxpy ( n-l+1, t, a(l,l), 1, a(l,j), 1 )
  520.  
  521.                 if ( (pl .le. j) .and. (j .le. pu) ) then
  522.  
  523.                   if ( qraux(j) /= 0.0D+00 ) then
  524.  
  525.                     tt = 1.0D+00 - ( abs ( a(l,j) ) / qraux(j) )**2
  526.                     tt = max ( tt, 0.0D+00 )
  527.                     t = tt
  528.                     tt = 1.0D+00 + 0.05D+00 * tt * (qraux(j)/work(j))**2
  529.  
  530.                     if ( tt /= 1.0D+00 ) then
  531.                       qraux(j) = qraux(j) * sqrt ( t )
  532.                     else
  533.                       qraux(j) = dnrm2 ( n-l, a(l+1,j), 1 )
  534.                       work(j) = qraux(j)
  535.                     end if
  536.  
  537.                   end if
  538.  
  539.                 end if
  540.  
  541.               end do
  542.  
  543. c  Save the transformation.
  544.  
  545.               qraux(l) = a(l,l)
  546.               a(l,l) = - nrmxl
  547.  
  548.             end if
  549.  
  550.           end if
  551.  
  552.         end do
  553.  
  554.         return
  555.       end
  556.       subroutine dqrls ( a, lda, m, n, tol, kr, b, x,r,jpvt,qraux,work,
  557.      #  itask, ind )
  558.  
  559. c*****************************************************************************80
  560. c
  561. c! DQRLS factors and solves a linear system in the least squares sense.
  562. c
  563. c  Discussion:
  564. c
  565. c    The linear system may be overdetermined, underdetermined or singular.
  566. c    The solution is obtained using a QR factorization of the
  567. c    coefficient matrix.
  568. c
  569. c    DQRLS can be efficiently used to solve several least squares
  570. c    problems with the same matrix A.  The first system is solved
  571. c    with ITASK = 1.  The subsequent systems are solved with
  572. c    ITASK = 2, to avoid the recomputation of the matrix factors.
  573. c    The parameters KR, JPVT, and QRAUX must not be modified
  574. c    between calls to DQRLS.
  575. c
  576. c    DQRLS is used to solve in a least squares sense
  577. c    overdetermined, underdetermined and singular linear systems.
  578. c    The system is A*X approximates B where A is M by N.
  579. c    B is a given M-vector, and X is the N-vector to be computed.
  580. c    A solution X is found which minimimzes the sum of squares (2-norm)
  581. c    of the residual,  A*X - B.
  582. c
  583. c    The numerical rank of A is determined using the tolerance TOL.
  584. c
  585. c    DQRLS uses the LINPACK subroutine DQRDC to compute the QR
  586. c    factorization, with column pivoting, of an M by N matrix A.
  587. c
  588. c  Modified:
  589. c
  590. c    15 April 2012
  591. c
  592. c  Reference:
  593. c
  594. c    David Kahaner, Cleve Moler, Steven Nash,
  595. c    Numerical Methods and Software,
  596. c    Prentice Hall, 1989,
  597. c    ISBN: 0-13-627258-4,
  598. c    LC: TA345.K34.
  599. c
  600. c  Parameters:
  601. c
  602. c    Input/output, real*8 A(LDA,N), an M by N matrix.
  603. c    On input, the matrix whose decomposition is to be computed.
  604. c    In a least squares data fitting problem, A(I,J) is the
  605. c    value of the J-th basis (model) function at the I-th data point.
  606. c    On output, A contains the output from DQRDC.  The triangular matrix R
  607. c    of the QR factorization is contained in the upper triangle and
  608. c    information needed to recover the orthogonal matrix Q is stored
  609. c    below the diagonal in A and in the vector QRAUX.
  610. c
  611. c    Input, integer LDA, the leading dimension of A.
  612. c
  613. c    Input, integer M, the number of rows of A.
  614. c
  615. c    Input, integer N, the number of columns of A.
  616. c
  617. c    Input, real*8 TOL, a relative tolerance used to determine the
  618. c    numerical rank.  The problem should be scaled so that all the elements
  619. c    of A have roughly the same absolute accuracy EPS.  Then a reasonable
  620. c    value for TOL is roughly EPS divided by the magnitude of the largest
  621. c    element.
  622. c
  623. c    Output, integer KR, the numerical rank.
  624. c
  625. c    Input, real*8 B(M), the right hand side of the linear system.
  626. c
  627. c    Output, real*8 X(N), a least squares solution to the linear
  628. c    system.
  629. c
  630. c    Output, real*8 R(M), the residual, B - A*X.  R may
  631. c    overwrite B.
  632. c
  633. c    Workspace, integer JPVT(N), required if ITASK = 1.
  634. c    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
  635. c    independent to within the tolerance TOL and the remaining columns
  636. c    are linearly dependent.  ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
  637. c    of the condition number of the matrix of independent columns,
  638. c    and of R.  This estimate will be .le. 1/TOL.
  639. c
  640. c    Workspace, real*8 QRAUX(N), required if ITASK = 1.
  641. c
  642. c    Workspace, real*8 WORK(N), required if ITASK = 1.
  643. c
  644. c    Input, integer ITASK.
  645. c    1, DQRLS factors the matrix A and solves the least squares problem.
  646. c    2, DQRLS assumes that the matrix A was factored with an earlier
  647. c       call to DQRLS, and only solves the least squares problem.
  648. c
  649. c    Output, integer IND, error code.
  650. c    0:  no error
  651. c    -1: LDA .lt. M   (fatal error)
  652. c    -2: N .lt. 1     (fatal error)
  653. c    -3: ITASK .lt. 1 (fatal error)
  654. c
  655.         implicit real*8 (a-h,o-z)
  656.         implicit integer (i-n)
  657.  
  658.         integer lda
  659.         integer m
  660.         integer n
  661.  
  662.         real*8 a(lda,n)
  663.         real*8 b(m)
  664.         integer ind
  665.         integer itask
  666.         integer jpvt(n)
  667.         integer kr
  668.         real*8 qraux(n)
  669.         real*8 r(m)
  670.         real*8 tol
  671.         real*8 work(n)
  672.         real*8 x(n)
  673.  
  674.         if ( lda .lt. m ) then
  675.           write ( *, '(a)' ) ' '
  676.           write ( *, '(a)' ) 'DQRLS - Fatal error!'
  677.           write ( *, '(a)' ) '  LDA .lt. M.'
  678.           stop
  679.         end if
  680.  
  681.         if ( n .le. 0 ) then
  682.           write ( *, '(a)' ) ' '
  683.           write ( *, '(a)' ) 'DQRLS - Fatal error!'
  684.           write ( *, '(a)' ) '  N .le. 0.'
  685.           stop
  686.         end if
  687.  
  688.         if ( itask .lt. 1 ) then
  689.           write ( *, '(a)' ) ' '
  690.           write ( *, '(a)' ) 'DQRLS - Fatal error!'
  691.           write ( *, '(a)' ) '  ITASK .lt. 1.'
  692.           stop
  693.         end if
  694.  
  695.         ind = 0
  696.  
  697. c  Factor the matrix.
  698.  
  699.         if ( itask .eq. 1 ) then
  700.           call dqrank ( a, lda, m, n, tol, kr, jpvt, qraux, work )
  701.         end if
  702.  
  703. c  Solve the least-squares problem.
  704.  
  705.         call dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
  706.  
  707.         return
  708.       end
  709.       subroutine dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
  710.  
  711. c*****************************************************************************80
  712. c
  713. c! DQRLSS solves a linear system in a least squares sense.
  714. c
  715. c  Discussion:
  716. c
  717. c    DQRLSS must be preceeded by a call to DQRANK.
  718. c
  719. c    The system is to be solved is
  720. c      A * X = B
  721. c    where
  722. c      A is an M by N matrix with rank KR, as determined by DQRANK,
  723. c      B is a given M-vector,
  724. c      X is the N-vector to be computed.
  725. c
  726. c    A solution X, with at most KR nonzero components, is found which
  727. c    minimizes the 2-norm of the residual (A*X-B).
  728. c
  729. c    Once the matrix A has been formed, DQRANK should be
  730. c    called once to decompose it.  Then, for each right hand
  731. c    side B, DQRLSS should be called once to obtain the
  732. c    solution and residual.
  733. c
  734. c  Modified:
  735. c
  736. c    15 April 2012
  737. c
  738. c  Parameters:
  739. c
  740. c    Input, real*8 A(LDA,N), the QR factorization information
  741. c    from DQRANK.  The triangular matrix R of the QR factorization is
  742. c    contained in the upper triangle and information needed to recover
  743. c    the orthogonal matrix Q is stored below the diagonal in A and in
  744. c    the vector QRAUX.
  745. c
  746. c    Input, integer LDA, the leading dimension of A, which must
  747. c    be at least M.
  748. c
  749. c    Input, integer M, the number of rows of A.
  750. c
  751. c    Input, integer N, the number of columns of A.
  752. c
  753. c    Input, integer KR, the rank of the matrix, as estimated
  754. c    by DQRANK.
  755. c
  756. c    Input, real*8 B(M), the right hand side of the linear system.
  757. c
  758. c    Output, real*8 X(N), a least squares solution to the
  759. c    linear system.
  760. c
  761. c    Output, real*8 R(M), the residual, B - A*X.  R may
  762. c    overwite B.
  763. c
  764. c    Input, integer JPVT(N), the pivot information from DQRANK.
  765. c    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
  766. c    independent to within the tolerance TOL and the remaining columns
  767. c    are linearly dependent.
  768. c
  769. c    Input, real*8 QRAUX(N), auxiliary information from DQRANK
  770. c    defining the QR factorization.
  771.  
  772.         implicit real*8 (a-h,o-z)
  773.         implicit integer (i-n)
  774.  
  775.         integer lda
  776.         integer m
  777.         integer n
  778.  
  779.         real*8 a(lda,n)
  780.         real*8 b(m)
  781.         integer info
  782.         integer j
  783.         integer job
  784.         integer jpvt(n)
  785.         integer k
  786.         integer kr
  787.         real*8 qraux(n)
  788.         real*8 r(m)
  789.         real*8 t
  790.         real*8 x(n)
  791.  
  792.         if ( kr /= 0 ) then
  793.           job = 110
  794.           call dqrsl ( a, lda, m, kr, qraux, b, r, r, x, r, r, job,info)
  795.         end if
  796.  
  797.         jpvt(1:n) = - jpvt(1:n)
  798.  
  799.         x(kr+1:n) = 0.0D+00
  800.  
  801.         do j = 1, n
  802.  
  803.           if ( jpvt(j) .le. 0 ) then
  804.  
  805.             k = -jpvt(j)
  806.             jpvt(j) = k
  807.  
  808.             do while ( k /= j )
  809.               t = x(j)
  810.               x(j) = x(k)
  811.               x(k) = t
  812.               jpvt(k) = -jpvt(k)
  813.               k = jpvt(k)
  814.             end do
  815.  
  816.           end if
  817.  
  818.         end do
  819.  
  820.         return
  821.       end
  822.       subroutine dqrsl ( a, lda, n, k,qraux,y,qy,qty,b,rsd,ab,job,info)
  823.  
  824. c*****************************************************************************80
  825. c
  826. c! DQRSL computes transformations, projections, and least squares solutions.
  827. c
  828. c  Discussion:
  829. c
  830. c    DQRSL requires the output of DQRDC.
  831. c
  832. c    For K .le. min(N,P), let AK be the matrix
  833. c
  834. c      AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
  835. c
  836. c    formed from columns JPVT(1), ..., JPVT(K) of the original
  837. c    N by P matrix A that was input to DQRDC.  If no pivoting was
  838. c    done, AK consists of the first K columns of A in their
  839. c    original order.  DQRDC produces a factored orthogonal matrix Q
  840. c    and an upper triangular matrix R such that
  841. c
  842. c      AK = Q * (R)
  843. c               (0)
  844. c
  845. c    This information is contained in coded form in the arrays
  846. c    A and QRAUX.
  847. c
  848. c    The parameters QY, QTY, B, RSD, and AB are not referenced
  849. c    if their computation is not requested and in this case
  850. c    can be replaced by dummy variables in the calling program.
  851. c    To save storage, the user may in some cases use the same
  852. c    array for different parameters in the calling sequence.  A
  853. c    frequently occuring example is when one wishes to compute
  854. c    any of B, RSD, or AB and does not need Y or QTY.  In this
  855. c    case one may identify Y, QTY, and one of B, RSD, or AB, while
  856. c    providing separate arrays for anything else that is to be
  857. c    computed.
  858. c
  859. c    Thus the calling sequence
  860. c
  861. c      call dqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
  862. c
  863. c    will result in the computation of B and RSD, with RSD
  864. c    overwriting Y.  More generally, each item in the following
  865. c    list contains groups of permissible identifications for
  866. c    a single calling sequence.
  867. c
  868. c      1. (Y,QTY,B) (RSD) (AB) (QY)
  869. c
  870. c      2. (Y,QTY,RSD) (B) (AB) (QY)
  871. c
  872. c      3. (Y,QTY,AB) (B) (RSD) (QY)
  873. c
  874. c      4. (Y,QY) (QTY,B) (RSD) (AB)
  875. c
  876. c      5. (Y,QY) (QTY,RSD) (B) (AB)
  877. c
  878. c      6. (Y,QY) (QTY,AB) (B) (RSD)
  879. c
  880. c    In any group the value returned in the array allocated to
  881. c    the group corresponds to the last member of the group.
  882. c
  883. c  Modified:
  884. c
  885. c    15 April 2012
  886. c
  887. c  Author:
  888. c
  889. c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart
  890. c
  891. c  Reference:
  892. c
  893. c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
  894. c    LINPACK User's Guide,
  895. c    SIAM, 1979,
  896. c    ISBN13: 978-0-898711-72-1,
  897. c    LC: QA214.L56.
  898. c
  899. c  Parameters:
  900. c
  901. c    Input, real*8 A(LDA,P), contains the output of DQRDC.
  902. c
  903. c    Input, integer LDA, the leading dimension of the array A.
  904. c
  905. c    Input, integer N, the number of rows of the matrix AK.  It 
  906. c    must have the same value as N in DQRDC.
  907. c
  908. c    Input, integer K, the number of columns of the matrix AK.  K
  909. c    must not be greater than min(N,P), where P is the same as in the
  910. c    calling sequence to DQRDC.
  911. c
  912. c    Input, real*8 QRAUX(P), the auxiliary output from DQRDC.
  913. c
  914. c    Input, real*8 Y(N), a vector to be manipulated by DQRSL.
  915. c
  916. c    Output, real*8 QY(N), contains Q * Y, if requested.
  917. c
  918. c    Output, real*8 QTY(N), contains Q' * Y, if requested.
  919. c
  920. c    Output, real*8 B(K), the solution of the least squares problem
  921. c      minimize norm2 ( Y - AK * B),
  922. c    if its computation has been requested.  Note that if pivoting was
  923. c    requested in DQRDC, the J-th component of B will be associated with
  924. c    column JPVT(J) of the original matrix A that was input into DQRDC.
  925. c
  926. c    Output, real*8 RSD(N), the least squares residual Y - AK * B,
  927. c    if its computation has been requested.  RSD is also the orthogonal
  928. c    projection of Y onto the orthogonal complement of the column space
  929. c    of AK.
  930. c
  931. c    Output, real*8 AB(N), the least squares approximation Ak * B,
  932. c    if its computation has been requested.  AB is also the orthogonal
  933. c    projection of Y onto the column space of A.
  934. c
  935. c    Input, integer JOB, specifies what is to be computed.  JOB has
  936. c    the decimal expansion ABCDE, with the following meaning:
  937. c
  938. c      if A /= 0, compute QY.
  939. c      if B /= 0, compute QTY.
  940. c      if C /= 0, compute QTY and B.
  941. c      if D /= 0, compute QTY and RSD.
  942. c      if E /= 0, compute QTY and AB.
  943. c
  944. c    Note that a request to compute B, RSD, or AB automatically triggers
  945. c    the computation of QTY, for which an array must be provided in the
  946. c    calling sequence.
  947. c
  948. c    Output, integer INFO, is zero unless the computation of B has
  949. c    been requested and R is exactly singular.  In this case, INFO is the
  950. c    index of the first zero diagonal element of R, and B is left unaltered.
  951.  
  952.         implicit real*8 (a-h,o-z)
  953.         implicit integer (i-n)
  954.  
  955.         integer k
  956.         integer lda
  957.         integer n
  958.  
  959.         real*8 a(lda,*)
  960.         real*8 ab(n)
  961.         real*8 b(k)
  962.         logical cab
  963.         logical cb
  964.         logical cqty
  965.         logical cqy
  966.         logical cr
  967.         integer info
  968.         integer j
  969.         integer jj
  970.         integer job
  971.         integer ju
  972.         integer kp1
  973.         real*8 qraux(*)
  974.         real*8 qty(n)
  975.         real*8 qy(n)
  976.         real*8 rsd(n)
  977.         real*8 ddot
  978.         real*8 t
  979.         real*8 temp
  980.         real*8 y(n)
  981.  
  982. c  set info flag.
  983.  
  984.         info = 0
  985.  
  986. c  Determine what is to be computed.
  987.  
  988.         cqy =        job / 10000         /= 0
  989.         cqty = mod ( job,  10000 )       /= 0
  990.         cb =   mod ( job,   1000 ) / 100 /= 0
  991.         cr =   mod ( job,    100 ) /  10 /= 0
  992.         cab =  mod ( job,     10 )       /= 0
  993.  
  994.         ju = min ( k, n-1 )
  995.  
  996. c  Special action when N = 1.
  997.  
  998.         if ( ju .eq. 0 ) then
  999.  
  1000.           if ( cqy ) then
  1001.             qy(1) = y(1)
  1002.           end if
  1003.  
  1004.           if ( cqty ) then
  1005.             qty(1) = y(1)
  1006.           end if
  1007.  
  1008.           if ( cab ) then
  1009.             ab(1) = y(1)
  1010.            end if
  1011.  
  1012.           if ( cb ) then
  1013.  
  1014.             if ( a(1,1) .eq. 0.0D+00 ) then
  1015.               info = 1
  1016.             else
  1017.               b(1) = y(1) / a(1,1)
  1018.             end if
  1019.  
  1020.           end if
  1021.  
  1022.           if ( cr ) then
  1023.             rsd(1) = 0.0D+00
  1024.           end if
  1025.  
  1026.           return
  1027.  
  1028.         end if
  1029.  
  1030. c  Set up to compute QY or QTY.
  1031.  
  1032.         if ( cqy ) then
  1033.           qy(1:n) = y(1:n)
  1034.         end if
  1035.  
  1036.         if ( cqty ) then
  1037.           qty(1:n) = y(1:n)
  1038.         end if
  1039.  
  1040. c  Compute QY.
  1041.  
  1042.         if ( cqy ) then
  1043.  
  1044.           do jj = 1, ju
  1045.  
  1046.             j = ju - jj + 1
  1047.  
  1048.             if ( qraux(j) /= 0.0D+00 ) then
  1049.               temp = a(j,j)
  1050.               a(j,j) = qraux(j)
  1051.               t = - ddot ( n-j+1, a(j,j), 1, qy(j), 1 ) / a(j,j)
  1052.               call daxpy ( n-j+1, t, a(j,j), 1, qy(j), 1 )
  1053.               a(j,j) = temp
  1054.             end if
  1055.  
  1056.           end do
  1057.  
  1058.         end if
  1059.  
  1060. c  Compute Q'*Y.
  1061.  
  1062.            if ( cqty ) then
  1063.  
  1064.               do j = 1, ju
  1065.                  if ( qraux(j) /= 0.0D+00 ) then
  1066.                     temp = a(j,j)
  1067.                     a(j,j) = qraux(j)
  1068.                     t = - ddot ( n-j+1, a(j,j), 1, qty(j), 1 ) / a(j,j)
  1069.                     call daxpy ( n-j+1, t, a(j,j), 1, qty(j), 1 )
  1070.                     a(j,j) = temp
  1071.                  end if
  1072.               end do
  1073.  
  1074.            end if
  1075.  
  1076. c  Set up to compute B, RSD, or AB.
  1077.  
  1078.            if ( cb ) then
  1079.              b(1:k) = qty(1:k)
  1080.            end if
  1081.  
  1082.            kp1 = k + 1
  1083.  
  1084.            if ( cab ) then
  1085.              ab(1:k) = qty(1:k)
  1086.            end if
  1087.  
  1088.            if ( cr .and. (k .lt. n) ) then
  1089.              rsd(k+1:n) = qty(k+1:n)
  1090.            end if
  1091.  
  1092.            if ( cab .and. (k+1 .le. n) ) then
  1093.               ab(k+1:n) = 0.0D+00
  1094.            end if
  1095.  
  1096.            if ( cr ) then
  1097.               rsd(1:k) = 0.0D+00
  1098.            end if
  1099.  
  1100. c  Compute B.
  1101.  
  1102.            if ( cb ) then
  1103.  
  1104.               do jj = 1, k
  1105.  
  1106.                  j = k - jj + 1
  1107.  
  1108.                  if ( a(j,j) .eq. 0.0D+00 ) then
  1109.                     info = j
  1110. c                    exit
  1111.                  go to 10
  1112.                  end if
  1113.  
  1114.                  b(j) = b(j)/a(j,j)
  1115.  
  1116.                  if ( j /= 1 ) then
  1117.                     t = -b(j)
  1118.                     call daxpy ( j-1, t, a(1,j), 1, b, 1 )
  1119.                  end if
  1120.  
  1121.               end do
  1122.    10      continue   
  1123.            end if
  1124.  
  1125.            if ( cr .or. cab ) then
  1126.  
  1127. c  Compute RSD or AB as required.
  1128.  
  1129.               do jj = 1, ju
  1130.  
  1131.                  j = ju - jj + 1
  1132.  
  1133.                  if ( qraux(j) /= 0.0D+00 ) then
  1134.  
  1135.                     temp = a(j,j)
  1136.                     a(j,j) = qraux(j)
  1137.  
  1138.                     if ( cr ) then
  1139.                      t = - ddot ( n-j+1, a(j,j), 1, rsd(j), 1 ) / a(j,j)
  1140.                      call daxpy ( n-j+1, t, a(j,j), 1, rsd(j), 1 )
  1141.                     end if
  1142.  
  1143.                     if ( cab ) then
  1144.                       t = - ddot ( n-j+1, a(j,j), 1, ab(j), 1 ) / a(j,j)
  1145.                        call daxpy ( n-j+1, t, a(j,j), 1, ab(j), 1 )
  1146.                     end if
  1147.  
  1148.                     a(j,j) = temp
  1149.  
  1150.                  end if
  1151.  
  1152.               end do
  1153.  
  1154.         end if
  1155.  
  1156.         return
  1157.       end
  1158.  
  1159.       subroutine dscal ( n, sa, x, incx )
  1160.  
  1161. c*****************************************************************************80
  1162. c
  1163. c! DSCAL scales a vector by a constant.
  1164. c
  1165. c  Discussion:
  1166. c
  1167. c    This routine uses double precision real arithmetic.
  1168. c
  1169. c  Licensing:
  1170. c
  1171. c    This code is distributed under the GNU LGPL license. 
  1172. c
  1173. c  Modified:
  1174. c
  1175. c    08 April 1999
  1176. c
  1177. c  Author:
  1178. c
  1179. c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
  1180. c    David Kincaid, Fred Krogh.
  1181. c    FORTRAN90 version by John Burkardt.
  1182. c
  1183. c  Reference:
  1184. c
  1185. c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
  1186. c    LINPACK User's Guide,
  1187. c    SIAM, 1979,
  1188. c    ISBN13: 978-0-898711-72-1,
  1189. c    LC: QA214.L56.
  1190. c
  1191. c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
  1192. c    Algorithm 539, 
  1193. c    Basic Linear Algebra Subprograms for Fortran Usage,
  1194. c    ACM Transactions on Mathematical Software,
  1195. c    Volume 5, Number 3, September 1979, pages 308-323.
  1196. c
  1197. c  Parameters:
  1198. c
  1199. c    Input, integer N, the number of entries in the vector.
  1200. c
  1201. c    Input, real*8 SA, the multiplier.
  1202. c
  1203. c    Input/output, real*8 X(*), the vector to be scaled.
  1204. c
  1205. c    Input, integer INCX, the increment between successive 
  1206. c    entries of X.
  1207. c
  1208.         implicit real*8 (a-h,o-z)
  1209.         implicit integer (i-n)
  1210.  
  1211.         integer i
  1212.         integer incx
  1213.         integer ix
  1214.         integer m
  1215.         integer n
  1216.         real*8 sa
  1217.         real*8 x(*)
  1218.  
  1219.         if ( n .le. 0 ) then
  1220.  
  1221.         else if ( incx .eq. 1 ) then
  1222.  
  1223.           m = mod ( n, 5 )
  1224.  
  1225.           x(1:m) = sa * x(1:m)
  1226.  
  1227.           do i = m+1, n, 5
  1228.             x(i)   = sa * x(i)
  1229.             x(i+1) = sa * x(i+1)
  1230.             x(i+2) = sa * x(i+2)
  1231.             x(i+3) = sa * x(i+3)
  1232.             x(i+4) = sa * x(i+4)
  1233.           end do
  1234.  
  1235.         else
  1236.  
  1237.           if ( 0 .le. incx ) then
  1238.             ix = 1
  1239.           else
  1240.             ix = ( - n + 1 ) * incx + 1
  1241.           end if
  1242.  
  1243.           do i = 1, n
  1244.             x(ix) = sa * x(ix)
  1245.             ix = ix + incx
  1246.           end do
  1247.  
  1248.         end if
  1249.  
  1250.         return
  1251.       end
  1252.  
  1253.       subroutine dswap ( n, x, incx, y, incy )
  1254.  
  1255. c*****************************************************************************80
  1256. c
  1257. c! DSWAP interchanges two vectors.
  1258. c
  1259. c  Discussion:
  1260. c
  1261. c    This routine uses double precision real arithmetic.
  1262. c
  1263. c  Licensing:
  1264. c
  1265. c    This code is distributed under the GNU LGPL license. 
  1266. c
  1267. c  Modified:
  1268. c
  1269. c    08 April 1999
  1270. c
  1271. c  Author:
  1272. c
  1273. c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
  1274. c    David Kincaid, Fred Krogh.
  1275. c    FORTRAN90 version by John Burkardt.
  1276. c
  1277. c  Reference:
  1278. c
  1279. c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
  1280. c    LINPACK User's Guide,
  1281. c    SIAM, 1979,
  1282. c    ISBN13: 978-0-898711-72-1,
  1283. c    LC: QA214.L56.
  1284. c
  1285. c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
  1286. c    Algorithm 539, 
  1287. c    Basic Linear Algebra Subprograms for Fortran Usage,
  1288. c    ACM Transactions on Mathematical Software, 
  1289. c    Volume 5, Number 3, September 1979, pages 308-323.
  1290. c
  1291. c  Parameters:
  1292. c
  1293. c    Input, integer N, the number of entries in the vectors.
  1294. c
  1295. c    Input/output, real*8 X(*), one of the vectors to swap.
  1296. c
  1297. c    Input, integer INCX, the increment between successive 
  1298. c    entries of X.
  1299. c
  1300. c    Input/output, real*8 Y(*), one of the vectors to swap.
  1301. c
  1302. c    Input, integer INCY, the increment between successive 
  1303. c    elements of Y.
  1304. c
  1305.         implicit real*8 (a-h,o-z)
  1306.         implicit integer (i-n)
  1307.  
  1308.         integer i
  1309.         integer incx
  1310.         integer incy
  1311.         integer ix
  1312.         integer iy
  1313.         integer m
  1314.         integer n
  1315.         real*8 temp
  1316.         real*8 x(*)
  1317.         real*8 y(*)
  1318.  
  1319.         if ( n .le. 0 ) then
  1320.  
  1321.         else if ( (incx .eq. 1) .and.( incy .eq. 1) ) then
  1322.  
  1323.           m = mod ( n, 3 )
  1324.  
  1325.           do i = 1, m
  1326.             temp = x(i)
  1327.             x(i) = y(i)
  1328.             y(i) = temp
  1329.           end do
  1330.  
  1331.           do i = m + 1, n, 3
  1332.  
  1333.             temp = x(i)
  1334.             x(i) = y(i)
  1335.             y(i) = temp
  1336.  
  1337.             temp = x(i+1)
  1338.             x(i+1) = y(i+1)
  1339.             y(i+1) = temp
  1340.  
  1341.             temp = x(i+2)
  1342.             x(i+2) = y(i+2)
  1343.             y(i+2) = temp
  1344.  
  1345.           end do
  1346.  
  1347.         else
  1348.  
  1349.           if ( 0 .le. incx ) then
  1350.             ix = 1
  1351.           else
  1352.             ix = ( - n + 1 ) * incx + 1
  1353.           end if
  1354.  
  1355.           if ( 0 .le. incy ) then
  1356.             iy = 1
  1357.           else
  1358.             iy = ( - n + 1 ) * incy + 1
  1359.           end if
  1360.  
  1361.           do i = 1, n
  1362.             temp = x(ix)
  1363.             x(ix) = y(iy)
  1364.             y(iy) = temp
  1365.             ix = ix + incx
  1366.             iy = iy + incy
  1367.           end do
  1368.  
  1369.         end if
  1370.  
  1371.         return
  1372.       end
  1373.  
  1374.  
  1375.  
  1376.  
  1377.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales