Télécharger la plaquette Cast3M


fr en

Télécharger dqrlss.eso

Retour à la liste

Numérotation des lignes : 

dqrlss

  1. C DQRLSS    SOURCE    PV090527  23/02/07    11:58:37     11592          
  2.       subroutine dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
  3.  
  4. c*****************************************************************************80
  5. c
  6. c! DQRLSS solves a linear system in a least squares sense.
  7. c
  8. c  Discussion:
  9. c
  10. c    DQRLSS must be preceeded by a call to DQRANK.
  11. c
  12. c    The system is to be solved is
  13. c      A * X = B
  14. c    where
  15. c      A is an M by N matrix with rank KR, as determined by DQRANK,
  16. c      B is a given M-vector,
  17. c      X is the N-vector to be computed.
  18. c
  19. c    A solution X, with at most KR nonzero components, is found which
  20. c    minimizes the 2-norm of the residual (A*X-B).
  21. c
  22. c    Once the matrix A has been formed, DQRANK should be
  23. c    called once to decompose it.  Then, for each right hand
  24. c    side B, DQRLSS should be called once to obtain the
  25. c    solution and residual.
  26. c
  27. c  Modified:
  28. c
  29. c    15 April 2012
  30. c
  31. c  Parameters:
  32. c
  33. c    Input, real*8 A(LDA,N), the QR factorization information
  34. c    from DQRANK.  The triangular matrix R of the QR factorization is
  35. c    contained in the upper triangle and information needed to recover
  36. c    the orthogonal matrix Q is stored below the diagonal in A and in
  37. c    the vector QRAUX.
  38. c
  39. c    Input, integer LDA, the leading dimension of A, which must
  40. c    be at least M.
  41. c
  42. c    Input, integer M, the number of rows of A.
  43. c
  44. c    Input, integer N, the number of columns of A.
  45. c
  46. c    Input, integer KR, the rank of the matrix, as estimated
  47. c    by DQRANK.
  48. c
  49. c    Input, real*8 B(M), the right hand side of the linear system.
  50. c
  51. c    Output, real*8 X(N), a least squares solution to the
  52. c    linear system.
  53. c
  54. c    Output, real*8 R(M), the residual, B - A*X.  R may
  55. c    overwite B.
  56. c
  57. c    Input, integer JPVT(N), the pivot information from DQRANK.
  58. c    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
  59. c    independent to within the tolerance TOL and the remaining columns
  60. c    are linearly dependent.
  61. c
  62. c    Input, real*8 QRAUX(N), auxiliary information from DQRANK
  63. c    defining the QR factorization.
  64.  
  65.         implicit real*8 (a-h,o-z)
  66.         implicit integer (i-n)
  67.  
  68.         integer lda
  69.         integer m
  70.         integer n
  71.  
  72.         real*8 a(lda,n)
  73.         real*8 b(m)
  74.         integer info
  75.         integer j
  76.         integer job
  77.         integer jpvt(n)
  78.         integer k
  79.         integer kr
  80.         real*8 qraux(n)
  81.         real*8 r(m)
  82.         real*8 t
  83.         real*8 x(n)
  84.  
  85.         if ( kr /= 0 ) then
  86.           job = 110
  87.           call dqrsl ( a, lda, m, kr, qraux, b, r, r, x, r, r, job,info)
  88.         end if
  89.  
  90.         jpvt(1:n) = - jpvt(1:n)
  91.  
  92.         x(kr+1:n) = 0.0D+00
  93.  
  94.         do j = 1, n
  95.  
  96.           if ( jpvt(j) .le. 0 ) then
  97.  
  98.             k = -jpvt(j)
  99.             jpvt(j) = k
  100.  
  101.             do while ( k /= j )
  102.               t = x(j)
  103.               x(j) = x(k)
  104.               x(k) = t
  105.               jpvt(k) = -jpvt(k)
  106.               k = jpvt(k)
  107.             end do
  108.  
  109.           end if
  110.  
  111.         end do
  112.  
  113.         return
  114.       end
  115.  
  116.  

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