home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / ludec / lapack.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  12.3 KB  |  420 lines
  1. /* *****************************************************************************
  2. *
  3. * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  4. * All Rights Reserved.
  5. *
  6. * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  7. * the contents of this file may not be disclosed to third parties, copied or
  8. * duplicated in any form, in whole or in part, without the prior written
  9. * permission of Silicon Graphics, Inc.
  10. *
  11. * RESTRICTED RIGHTS LEGEND:
  12. * Use, duplication or disclosure by the Government is subject to restrictions
  13. * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  14. * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  15. * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  16. * rights reserved under the Copyright Laws of the United States.
  17. *
  18. ***************************************************************************** */
  19.       SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
  20. *
  21. *  -- LAPACK routine (preliminary version) --
  22. *     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
  23. *     Courant Institute, NAG Ltd., and Rice University
  24. *     March 26, 1990
  25. *
  26. *     .. Scalar Arguments ..
  27.       INTEGER            INFO, LDA, M, N
  28. *     ..
  29. *     .. Array Arguments ..
  30.       INTEGER            IPIV( * )
  31.       DOUBLE PRECISION   A( LDA, * )
  32. *     ..
  33. *
  34. *  Purpose
  35. *  =======
  36. *
  37. *  DGETRF computes an LU factorization of a general m-by-n matrix A
  38. *  using partial pivoting with row interchanges.
  39. *
  40. *  The factorization has the form
  41. *     A = P * L * U
  42. *  where P is a permutation matrix, L is lower triangular with unit
  43. *  diagonal elements (lower trapezoidal if m > n), and U is upper
  44. *  triangular (upper trapezoidal if m < n).
  45. *
  46. *  This is the Level 3 BLAS version of the Crout algorithm.
  47. *
  48. *  Arguments
  49. *  =========
  50. *
  51. *  M       (input) INTEGER
  52. *          The number of rows of the matrix A.  M >= 0.
  53. *
  54. *  N       (input) INTEGER
  55. *          The number of columns of the matrix A.  N >= 0.
  56. *
  57. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  58. *          On entry, the m x n matrix to be factored.
  59. *          On exit, the factors L and U; the unit diagonal elements of L
  60. *          are not stored.
  61. *
  62. *  LDA     (input) INTEGER
  63. *          The leading dimension of the array A.  LDA >= max(1,M).
  64. *
  65. *  IPIV    (output) INTEGER array, dimension (min(M,N))
  66. *          The pivot indices.  Row i of the matrix was interchanged with
  67. *          row IPIV(i).
  68. *
  69. *  INFO    (output) INTEGER
  70. *          = 0: successful exit
  71. *          < 0: if INFO = -k, the k-th argument had an illegal value
  72. *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
  73. *               has been completed, but the factor U is exactly
  74. *               singular, and division by zero will occur if it is used
  75. *               to solve a system of equations or to compute the inverse
  76. *               of A.
  77. *
  78. *  =====================================================================
  79. *
  80. *     .. Parameters ..
  81.       DOUBLE PRECISION   ONE
  82.       PARAMETER          ( ONE = 1.0D+0 )
  83. *     ..
  84. *     .. Local Scalars ..
  85.       INTEGER            I, J, JB, NB
  86. *     ..
  87. *     .. External Subroutines ..
  88. c     EXTERNAL           DGEMM, DGETF2, DLASWP, DTRSM, ENVIR, XERBLA
  89.       EXTERNAL           DGEMM, DGETF2, DLASWP, DTRSM, XERBLA
  90. *     ..
  91. *     .. Intrinsic Functions ..
  92.       INTRINSIC          MAX, MIN
  93.  
  94. c     ..jpp -> Blocksize is fixed to 32 for demo
  95.       nb = 32
  96. *     ..
  97. *     .. Executable Statements ..
  98. *
  99. *     Test the input parameters.
  100. *
  101.       INFO = 0
  102.       IF( M.LT.0 ) THEN
  103.          INFO = -1
  104.       ELSE IF( N.LT.0 ) THEN
  105.          INFO = -2
  106.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  107.          INFO = -4
  108.       END IF
  109.       IF( INFO.NE.0 ) THEN
  110.          CALL XERBLA( 'DGETRF', -INFO )
  111.          RETURN
  112.       END IF
  113. *
  114. *     Quick return if possible.
  115. *
  116.       IF( M.EQ.0 .OR. N.EQ.0 )
  117.      $   RETURN
  118. *
  119. *     Determine the block size for this environment.
  120. *
  121. c     CALL ENVIR( 'Block', NB )
  122.       IF( NB.EQ.1 ) THEN
  123.          CALL DGETF2( M, N, A, LDA, IPIV, INFO )
  124.       ELSE
  125. *
  126.          DO 20 J = 1, MIN( M, N ), NB
  127.         call update(j-1)
  128.             JB = MIN( MIN( M, N )-J+1, NB )
  129. *
  130. *           Update diagonal and subdiagonal blocks.
  131. *
  132.             CALL DGEMM( 'No transpose', 'No transpose', M-J+1, JB, J-1,
  133.      $                  -ONE, A( J, 1 ), LDA, A( 1, J ), LDA, ONE,
  134.      $                  A( J, J ), LDA )
  135. *
  136. *           Factorize diagonal and subdiagonal blocks and test for
  137. *           singularity.
  138. *
  139.             CALL DGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), INFO )
  140. *
  141. *           Update pivot indices and apply the interchanges to the
  142. *           columns on either side of the current block.
  143. *
  144.             DO 10 I = J, MIN( M, J+JB-1 )
  145.                IPIV( I ) = J - 1 + IPIV( I )
  146.    10       CONTINUE
  147.             CALL DLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
  148.             CALL DLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1, IPIV,
  149.      $                   1 )
  150. *
  151.             IF( INFO.NE.0 )
  152.      $         GO TO 30
  153. *
  154.             IF( J+JB.LE.N ) THEN
  155. *
  156. *              Compute block row of U.
  157. *
  158.                CALL DGEMM( 'No transpose', 'No transpose', JB, N-J-JB+1,
  159.      $                     J-1, -ONE, A( J, 1 ), LDA, A( 1, J+JB ), LDA,
  160.      $                     ONE, A( J, J+JB ), LDA )
  161.                CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
  162.      $                     N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
  163.      $                     LDA )
  164.             END IF
  165.    20    CONTINUE
  166.     call update(m)
  167.       END IF
  168.       RETURN
  169. *
  170.    30 CONTINUE
  171.       INFO = INFO + J - 1
  172.       RETURN
  173. *
  174. *     End of DGETRF
  175. *
  176.       END
  177.       SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
  178. *
  179. *  -- LAPACK routine (preliminary version) --
  180. *     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
  181. *     Courant Institute, NAG Ltd., and Rice University
  182. *     March 26, 1990
  183. *
  184. *     .. Scalar Arguments ..
  185.       INTEGER            INFO, LDA, M, N
  186. *     ..
  187. *     .. Array Arguments ..
  188.       INTEGER            IPIV( * )
  189.       DOUBLE PRECISION   A( LDA, * )
  190. *     ..
  191. *
  192. *  Purpose
  193. *  =======
  194. *
  195. *  DGETF2 computes an LU factorization of a general m-by-n matrix A
  196. *  using partial pivoting with row interchanges.
  197. *
  198. *  The factorization has the form
  199. *     A = P * L * U
  200. *  where P is a permutation matrix, L is lower triangular with unit
  201. *  diagonal elements (lower trapezoidal if m > n), and U is upper
  202. *  triangular (upper trapezoidal if m < n).
  203. *
  204. *  This is the Level 2 BLAS version of the Crout algorithm.
  205. *
  206. *  Arguments
  207. *  =========
  208. *
  209. *  M       (input) INTEGER
  210. *          The number of rows of the matrix A.  M >= 0.
  211. *
  212. *  N       (input) INTEGER
  213. *          The number of columns of the matrix A.  N >= 0.
  214. *
  215. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  216. *          On entry, the m x n matrix to be factored.
  217. *          On exit, the factors L and U; the unit diagonal elements of L
  218. *          are not stored.
  219. *
  220. *  LDA     (input) INTEGER
  221. *          The leading dimension of the array A.  LDA >= max(1,M).
  222. *
  223. *  IPIV    (output) INTEGER array, dimension (min(M,N))
  224. *          The pivot indices.  Row i of the matrix was interchanged with
  225. *          row IPIV(i).
  226. *
  227. *  INFO    (output) INTEGER
  228. *          = 0: successful exit
  229. *          < 0: if INFO = -k, the k-th argument had an illegal value
  230. *          > 0: if INFO = k, U(k,k) is exactly zero. The factorization
  231. *               has been completed, but the factor U is exactly
  232. *               singular, and division by zero will occur if it is used
  233. *               to solve a system of equations or to compute the inverse
  234. *               of A.
  235. *
  236. *  =====================================================================
  237. *
  238. *     .. Parameters ..
  239.       DOUBLE PRECISION   ONE, ZERO
  240.       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  241. *     ..
  242. *     .. Local Scalars ..
  243.       INTEGER            J, JP
  244. *     ..
  245. *     .. External Functions ..
  246.       INTEGER            IDAMAX
  247.       EXTERNAL           IDAMAX
  248. *     ..
  249. *     .. External Subroutines ..
  250.       EXTERNAL           DGEMV, DSCAL, DSWAP, XERBLA
  251. *     ..
  252. *     .. Intrinsic Functions ..
  253.       INTRINSIC          MAX, MIN
  254. *     ..
  255. *     .. Executable Statements ..
  256. *
  257. *
  258. *     Test the input parameters.
  259. *
  260.       INFO = 0
  261.       IF( M.LT.0 ) THEN
  262.          INFO = -1
  263.       ELSE IF( N.LT.0 ) THEN
  264.          INFO = -2
  265.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  266.          INFO = -4
  267.       END IF
  268.       IF( INFO.NE.0 ) THEN
  269.          CALL XERBLA( 'DGETF2', -INFO )
  270.          RETURN
  271.       END IF
  272. *
  273. *     Quick return if possible.
  274. *
  275.       IF( M.EQ.0 .OR. N.EQ.0 )
  276.      $   RETURN
  277. *
  278.       DO 10 J = 1, MIN( M, N )
  279. *
  280. *        Update diagonal and subdiagonal elements in column J.
  281. *
  282.          CALL DGEMV( 'No transpose', M-J+1, J-1, -ONE, A( J, 1 ), LDA,
  283.      $               A( 1, J ), 1, ONE, A( J, J ), 1 )
  284. *
  285. *        Find pivot and test for singularity.
  286. *
  287.          JP = J - 1 + IDAMAX( M-J+1, A( J, J ), 1 )
  288.          IPIV( J ) = JP
  289.          IF( A( JP, J ).NE.ZERO ) THEN
  290. *
  291. *           Apply interchange to columns 1:J.
  292. *
  293.             IF( JP.NE.J )
  294.      $         CALL DSWAP( N, A( J, 1 ), LDA, A( JP, 1 ), LDA )
  295. *
  296. *           Compute elements J+1:M of J-th column.
  297. *
  298.             IF( J.LT.M )
  299.      $         CALL DSCAL( M-J, ONE/A( J, J ), A( J+1, J ), 1 )
  300.          ELSE
  301. *
  302. *           If A( JP, J ) is zero, set INFO to indicate that a zero
  303. *           pivot has been found.  INFO returns the index of the
  304. *           last zero pivot to the calling routine if there is more
  305. *           than one.
  306. *
  307.             INFO = J
  308.          END IF
  309. *
  310.          IF( J+1.LE.N ) THEN
  311. *
  312. *           Compute block row of U.
  313. *
  314.             CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
  315.      $                  A( J, 1 ), LDA, ONE, A( J, J+1 ), LDA )
  316.          END IF
  317.    10 CONTINUE
  318.       RETURN
  319. *
  320. *     End of DGETF2
  321. *
  322.       END
  323.       SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX )
  324. *
  325. *  -- LAPACK auxiliary routine (preliminary version) --
  326. *     Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
  327. *     Courant Institute, NAG Ltd., and Rice University
  328. *     March 26, 1990
  329. *
  330. *     .. Scalar Arguments ..
  331.       INTEGER            INCX, K1, K2, LDA, N
  332. *     ..
  333. *     .. Array Arguments ..
  334.       INTEGER            IPIV( * )
  335.       DOUBLE PRECISION   A( LDA, * )
  336. *     ..
  337. *
  338. *  Purpose
  339. *  =======
  340. *
  341. *  DLASWP performs a series of row interchanges on the M x N matrix A.
  342. *  One row interchange is initiated for each of rows K1 through K2 of A.
  343. *
  344. *  Arguments
  345. *  =========
  346. *
  347. *  N       (input) INTEGER
  348. *          The number of columns of the matrix A.
  349. *
  350. *  A       (input/output) DOUBLE PRECISION array, dimension( LDA, N )
  351. *          On entry, the M x N matrix to which the row interchanges will
  352. *          be applied.
  353. *          On exit, the permuted matrix.
  354. *
  355. *  LDA     (input) INTEGER
  356. *          The leading dimension of the array A.
  357. *
  358. *  K1      (input) INTEGER
  359. *          The first element of IPIV for which a row interchange will
  360. *          be done.
  361. *
  362. *  K2      (input) INTEGER
  363. *          The last element of IPIV for which a row interchange will
  364. *          be done.
  365. *
  366. *  IPIV    (input) INTEGER array, dimension( M*abs(INCX) )
  367. *          The vector of pivot indices.  Only the elements in positions
  368. *          K1 through K2 of IPIV are accessed.
  369. *          IPIV(K) = L implies rows K and L are to be interchanged.
  370. *
  371. *  INCX    (input) INTEGER
  372. *          The increment between succesive values of IPIV.  If IPIV
  373. *          is negative, the pivots are applied in reverse order.
  374. *
  375. *
  376. *     .. Local Scalars ..
  377.       INTEGER            I, IP, IX
  378. *     ..
  379. *     .. External Subroutines ..
  380.       EXTERNAL           DSWAP
  381. *     ..
  382. *     .. Executable Statements ..
  383. *
  384. *     Interchange row I with row IPIV(I) for each of rows K1 through K2.
  385. *
  386.       IF( INCX.EQ.0 )
  387.      $   RETURN
  388.       IF( INCX.GT.0 ) THEN
  389.          IX = K1
  390.       ELSE
  391.          IX = 1 + ( 1-K2 )*INCX
  392.       END IF
  393.       IF( INCX.EQ.1 ) THEN
  394.          DO 10 I = K1, K2
  395.             IP = IPIV( I )
  396.             IF( IP.NE.I )
  397.      $         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
  398.    10    CONTINUE
  399.       ELSE IF( INCX.GT.1 ) THEN
  400.          DO 20 I = K1, K2
  401.             IP = IPIV( IX )
  402.             IF( IP.NE.I )
  403.      $         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
  404.             IX = IX + INCX
  405.    20    CONTINUE
  406.       ELSE IF( INCX.LT.0 ) THEN
  407.          DO 30 I = K2, K1, -1
  408.             IP = IPIV( IX )
  409.             IF( IP.NE.I )
  410.      $         CALL DSWAP( N, A( I, 1 ), LDA, A( IP, 1 ), LDA )
  411.             IX = IX + INCX
  412.    30    CONTINUE
  413.       END IF
  414. *
  415.       RETURN
  416. *
  417. *     End of DLASWP
  418. *
  419.       END
  420.