LAPACK: TESTING/LIN/zlavsp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLAVSP( UPLO, TRANS, DIAG, N, NRHS, A, IPIV, B, LDB,
00002      $                   INFO )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          DIAG, TRANS, UPLO
00010       INTEGER            INFO, LDB, N, NRHS
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IPIV( * )
00014       COMPLEX*16         A( * ), B( LDB, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *     ZLAVSP  performs one of the matrix-vector operations
00021 *        x := A*x  or  x := A^T*x,
00022 *     where x is an N element vector and  A is one of the factors
00023 *     from the symmetric factorization computed by ZSPTRF.
00024 *     ZSPTRF produces a factorization of the form
00025 *          U * D * U^T     or     L * D * L^T,
00026 *     where U (or L) is a product of permutation and unit upper (lower)
00027 *     triangular matrices, U^T (or L^T) is the transpose of
00028 *     U (or L), and D is symmetric and block diagonal with 1 x 1 and
00029 *     2 x 2 diagonal blocks.  The multipliers for the transformations
00030 *     and the upper or lower triangular parts of the diagonal blocks
00031 *     are stored columnwise in packed format in the linear array A.
00032 *
00033 *     If TRANS = 'N' or 'n', ZLAVSP multiplies either by U or U * D
00034 *     (or L or L * D).
00035 *     If TRANS = 'C' or 'c', ZLAVSP multiplies either by U^T or D * U^T
00036 *     (or L^T or D * L^T ).
00037 *
00038 *  Arguments
00039 *  ==========
00040 *
00041 *  UPLO   - CHARACTER*1
00042 *           On entry, UPLO specifies whether the triangular matrix
00043 *           stored in A is upper or lower triangular.
00044 *              UPLO = 'U' or 'u'   The matrix is upper triangular.
00045 *              UPLO = 'L' or 'l'   The matrix is lower triangular.
00046 *           Unchanged on exit.
00047 *
00048 *  TRANS  - CHARACTER*1
00049 *           On entry, TRANS specifies the operation to be performed as
00050 *           follows:
00051 *              TRANS = 'N' or 'n'   x := A*x.
00052 *              TRANS = 'T' or 't'   x := A^T*x.
00053 *           Unchanged on exit.
00054 *
00055 *  DIAG   - CHARACTER*1
00056 *           On entry, DIAG specifies whether the diagonal blocks are
00057 *           assumed to be unit matrices, as follows:
00058 *              DIAG = 'U' or 'u'   Diagonal blocks are unit matrices.
00059 *              DIAG = 'N' or 'n'   Diagonal blocks are non-unit.
00060 *           Unchanged on exit.
00061 *
00062 *  N      - INTEGER
00063 *           On entry, N specifies the order of the matrix A.
00064 *           N must be at least zero.
00065 *           Unchanged on exit.
00066 *
00067 *  NRHS   - INTEGER
00068 *           On entry, NRHS specifies the number of right hand sides,
00069 *           i.e., the number of vectors x to be multiplied by A.
00070 *           NRHS must be at least zero.
00071 *           Unchanged on exit.
00072 *
00073 *  A      - COMPLEX*16 array, dimension( N*(N+1)/2 )
00074 *           On entry, A contains a block diagonal matrix and the
00075 *           multipliers of the transformations used to obtain it,
00076 *           stored as a packed triangular matrix.
00077 *           Unchanged on exit.
00078 *
00079 *  IPIV   - INTEGER array, dimension( N )
00080 *           On entry, IPIV contains the vector of pivot indices as
00081 *           determined by ZSPTRF.
00082 *           If IPIV( K ) = K, no interchange was done.
00083 *           If IPIV( K ) <> K but IPIV( K ) > 0, then row K was inter-
00084 *           changed with row IPIV( K ) and a 1 x 1 pivot block was used.
00085 *           If IPIV( K ) < 0 and UPLO = 'U', then row K-1 was exchanged
00086 *           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
00087 *           If IPIV( K ) < 0 and UPLO = 'L', then row K+1 was exchanged
00088 *           with row | IPIV( K ) | and a 2 x 2 pivot block was used.
00089 *
00090 *  B      - COMPLEX*16 array, dimension( LDB, NRHS )
00091 *           On entry, B contains NRHS vectors of length N.
00092 *           On exit, B is overwritten with the product A * B.
00093 *
00094 *  LDB    - INTEGER
00095 *           On entry, LDB contains the leading dimension of B as
00096 *           declared in the calling program.  LDB must be at least
00097 *           max( 1, N ).
00098 *           Unchanged on exit.
00099 *
00100 *  INFO   - INTEGER
00101 *           INFO is the error flag.
00102 *           On exit, a value of 0 indicates a successful exit.
00103 *           A negative value, say -K, indicates that the K-th argument
00104 *           has an illegal value.
00105 *
00106 *  =====================================================================
00107 *
00108 *     .. Parameters ..
00109       COMPLEX*16         ONE
00110       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
00111 *     ..
00112 *     .. Local Scalars ..
00113       LOGICAL            NOUNIT
00114       INTEGER            J, K, KC, KCNEXT, KP
00115       COMPLEX*16         D11, D12, D21, D22, T1, T2
00116 *     ..
00117 *     .. External Functions ..
00118       LOGICAL            LSAME
00119       EXTERNAL           LSAME
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           XERBLA, ZGEMV, ZGERU, ZSCAL, ZSWAP
00123 *     ..
00124 *     .. Intrinsic Functions ..
00125       INTRINSIC          ABS, MAX
00126 *     ..
00127 *     .. Executable Statements ..
00128 *
00129 *     Test the input parameters.
00130 *
00131       INFO = 0
00132       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00133          INFO = -1
00134       ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
00135      $          THEN
00136          INFO = -2
00137       ELSE IF( .NOT.LSAME( DIAG, 'U' ) .AND. .NOT.LSAME( DIAG, 'N' ) )
00138      $          THEN
00139          INFO = -3
00140       ELSE IF( N.LT.0 ) THEN
00141          INFO = -4
00142       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00143          INFO = -8
00144       END IF
00145       IF( INFO.NE.0 ) THEN
00146          CALL XERBLA( 'ZLAVSP ', -INFO )
00147          RETURN
00148       END IF
00149 *
00150 *     Quick return if possible.
00151 *
00152       IF( N.EQ.0 )
00153      $   RETURN
00154 *
00155       NOUNIT = LSAME( DIAG, 'N' )
00156 *------------------------------------------
00157 *
00158 *     Compute  B := A * B  (No transpose)
00159 *
00160 *------------------------------------------
00161       IF( LSAME( TRANS, 'N' ) ) THEN
00162 *
00163 *        Compute  B := U*B
00164 *        where U = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00165 *
00166          IF( LSAME( UPLO, 'U' ) ) THEN
00167 *
00168 *        Loop forward applying the transformations.
00169 *
00170             K = 1
00171             KC = 1
00172    10       CONTINUE
00173             IF( K.GT.N )
00174      $         GO TO 30
00175 *
00176 *           1 x 1 pivot block
00177 *
00178             IF( IPIV( K ).GT.0 ) THEN
00179 *
00180 *              Multiply by the diagonal element if forming U * D.
00181 *
00182                IF( NOUNIT )
00183      $            CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00184 *
00185 *              Multiply by P(K) * inv(U(K))  if K > 1.
00186 *
00187                IF( K.GT.1 ) THEN
00188 *
00189 *                 Apply the transformation.
00190 *
00191                   CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
00192      $                        LDB, B( 1, 1 ), LDB )
00193 *
00194 *                 Interchange if P(K) != I.
00195 *
00196                   KP = IPIV( K )
00197                   IF( KP.NE.K )
00198      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00199                END IF
00200                KC = KC + K
00201                K = K + 1
00202             ELSE
00203 *
00204 *              2 x 2 pivot block
00205 *
00206                KCNEXT = KC + K
00207 *
00208 *              Multiply by the diagonal block if forming U * D.
00209 *
00210                IF( NOUNIT ) THEN
00211                   D11 = A( KCNEXT-1 )
00212                   D22 = A( KCNEXT+K )
00213                   D12 = A( KCNEXT+K-1 )
00214                   D21 = D12
00215                   DO 20 J = 1, NRHS
00216                      T1 = B( K, J )
00217                      T2 = B( K+1, J )
00218                      B( K, J ) = D11*T1 + D12*T2
00219                      B( K+1, J ) = D21*T1 + D22*T2
00220    20             CONTINUE
00221                END IF
00222 *
00223 *              Multiply by  P(K) * inv(U(K))  if K > 1.
00224 *
00225                IF( K.GT.1 ) THEN
00226 *
00227 *                 Apply the transformations.
00228 *
00229                   CALL ZGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
00230      $                        LDB, B( 1, 1 ), LDB )
00231                   CALL ZGERU( K-1, NRHS, ONE, A( KCNEXT ), 1,
00232      $                        B( K+1, 1 ), LDB, B( 1, 1 ), LDB )
00233 *
00234 *                 Interchange if P(K) != I.
00235 *
00236                   KP = ABS( IPIV( K ) )
00237                   IF( KP.NE.K )
00238      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00239                END IF
00240                KC = KCNEXT + K + 1
00241                K = K + 2
00242             END IF
00243             GO TO 10
00244    30       CONTINUE
00245 *
00246 *        Compute  B := L*B
00247 *        where L = P(1)*inv(L(1))* ... *P(m)*inv(L(m)) .
00248 *
00249          ELSE
00250 *
00251 *           Loop backward applying the transformations to B.
00252 *
00253             K = N
00254             KC = N*( N+1 ) / 2 + 1
00255    40       CONTINUE
00256             IF( K.LT.1 )
00257      $         GO TO 60
00258             KC = KC - ( N-K+1 )
00259 *
00260 *           Test the pivot index.  If greater than zero, a 1 x 1
00261 *           pivot was used, otherwise a 2 x 2 pivot was used.
00262 *
00263             IF( IPIV( K ).GT.0 ) THEN
00264 *
00265 *              1 x 1 pivot block:
00266 *
00267 *              Multiply by the diagonal element if forming L * D.
00268 *
00269                IF( NOUNIT )
00270      $            CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00271 *
00272 *              Multiply by  P(K) * inv(L(K))  if K < N.
00273 *
00274                IF( K.NE.N ) THEN
00275                   KP = IPIV( K )
00276 *
00277 *                 Apply the transformation.
00278 *
00279                   CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00280      $                        LDB, B( K+1, 1 ), LDB )
00281 *
00282 *                 Interchange if a permutation was applied at the
00283 *                 K-th step of the factorization.
00284 *
00285                   IF( KP.NE.K )
00286      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00287                END IF
00288                K = K - 1
00289 *
00290             ELSE
00291 *
00292 *              2 x 2 pivot block:
00293 *
00294                KCNEXT = KC - ( N-K+2 )
00295 *
00296 *              Multiply by the diagonal block if forming L * D.
00297 *
00298                IF( NOUNIT ) THEN
00299                   D11 = A( KCNEXT )
00300                   D22 = A( KC )
00301                   D21 = A( KCNEXT+1 )
00302                   D12 = D21
00303                   DO 50 J = 1, NRHS
00304                      T1 = B( K-1, J )
00305                      T2 = B( K, J )
00306                      B( K-1, J ) = D11*T1 + D12*T2
00307                      B( K, J ) = D21*T1 + D22*T2
00308    50             CONTINUE
00309                END IF
00310 *
00311 *              Multiply by  P(K) * inv(L(K))  if K < N.
00312 *
00313                IF( K.NE.N ) THEN
00314 *
00315 *                 Apply the transformation.
00316 *
00317                   CALL ZGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00318      $                        LDB, B( K+1, 1 ), LDB )
00319                   CALL ZGERU( N-K, NRHS, ONE, A( KCNEXT+2 ), 1,
00320      $                        B( K-1, 1 ), LDB, B( K+1, 1 ), LDB )
00321 *
00322 *                 Interchange if a permutation was applied at the
00323 *                 K-th step of the factorization.
00324 *
00325                   KP = ABS( IPIV( K ) )
00326                   IF( KP.NE.K )
00327      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00328                END IF
00329                KC = KCNEXT
00330                K = K - 2
00331             END IF
00332             GO TO 40
00333    60       CONTINUE
00334          END IF
00335 *-------------------------------------------------
00336 *
00337 *     Compute  B := A^T * B  (transpose)
00338 *
00339 *-------------------------------------------------
00340       ELSE
00341 *
00342 *        Form  B := U^T*B
00343 *        where U  = P(m)*inv(U(m))* ... *P(1)*inv(U(1))
00344 *        and   U^T = inv(U^T(1))*P(1)* ... *inv(U^T(m))*P(m)
00345 *
00346          IF( LSAME( UPLO, 'U' ) ) THEN
00347 *
00348 *           Loop backward applying the transformations.
00349 *
00350             K = N
00351             KC = N*( N+1 ) / 2 + 1
00352    70       CONTINUE
00353             IF( K.LT.1 )
00354      $         GO TO 90
00355             KC = KC - K
00356 *
00357 *           1 x 1 pivot block.
00358 *
00359             IF( IPIV( K ).GT.0 ) THEN
00360                IF( K.GT.1 ) THEN
00361 *
00362 *                 Interchange if P(K) != I.
00363 *
00364                   KP = IPIV( K )
00365                   IF( KP.NE.K )
00366      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00367 *
00368 *                 Apply the transformation:
00369 *                    y := y - B' * conjg(x)
00370 *                 where x is a column of A and y is a row of B.
00371 *
00372                   CALL ZGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
00373      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00374                END IF
00375                IF( NOUNIT )
00376      $            CALL ZSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00377                K = K - 1
00378 *
00379 *           2 x 2 pivot block.
00380 *
00381             ELSE
00382                KCNEXT = KC - ( K-1 )
00383                IF( K.GT.2 ) THEN
00384 *
00385 *                 Interchange if P(K) != I.
00386 *
00387                   KP = ABS( IPIV( K ) )
00388                   IF( KP.NE.K-1 )
00389      $               CALL ZSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00390      $                           LDB )
00391 *
00392 *                 Apply the transformations.
00393 *
00394                   CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00395      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00396 *
00397                   CALL ZGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00398      $                        A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
00399                END IF
00400 *
00401 *              Multiply by the diagonal block if non-unit.
00402 *
00403                IF( NOUNIT ) THEN
00404                   D11 = A( KC-1 )
00405                   D22 = A( KC+K-1 )
00406                   D12 = A( KC+K-2 )
00407                   D21 = D12
00408                   DO 80 J = 1, NRHS
00409                      T1 = B( K-1, J )
00410                      T2 = B( K, J )
00411                      B( K-1, J ) = D11*T1 + D12*T2
00412                      B( K, J ) = D21*T1 + D22*T2
00413    80             CONTINUE
00414                END IF
00415                KC = KCNEXT
00416                K = K - 2
00417             END IF
00418             GO TO 70
00419    90       CONTINUE
00420 *
00421 *        Form  B := L^T*B
00422 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00423 *        and   L^T = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
00424 *
00425          ELSE
00426 *
00427 *           Loop forward applying the L-transformations.
00428 *
00429             K = 1
00430             KC = 1
00431   100       CONTINUE
00432             IF( K.GT.N )
00433      $         GO TO 120
00434 *
00435 *           1 x 1 pivot block
00436 *
00437             IF( IPIV( K ).GT.0 ) THEN
00438                IF( K.LT.N ) THEN
00439 *
00440 *                 Interchange if P(K) != I.
00441 *
00442                   KP = IPIV( K )
00443                   IF( KP.NE.K )
00444      $               CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00445 *
00446 *                 Apply the transformation
00447 *
00448                   CALL ZGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
00449      $                        LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00450                END IF
00451                IF( NOUNIT )
00452      $            CALL ZSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00453                KC = KC + N - K + 1
00454                K = K + 1
00455 *
00456 *           2 x 2 pivot block.
00457 *
00458             ELSE
00459                KCNEXT = KC + N - K + 1
00460                IF( K.LT.N-1 ) THEN
00461 *
00462 *              Interchange if P(K) != I.
00463 *
00464                   KP = ABS( IPIV( K ) )
00465                   IF( KP.NE.K+1 )
00466      $               CALL ZSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00467      $                           LDB )
00468 *
00469 *                 Apply the transformation
00470 *
00471                   CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
00472      $                        B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
00473      $                        B( K+1, 1 ), LDB )
00474 *
00475                   CALL ZGEMV( 'Transpose', N-K-1, NRHS, ONE,
00476      $                        B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
00477      $                        B( K, 1 ), LDB )
00478                END IF
00479 *
00480 *              Multiply by the diagonal block if non-unit.
00481 *
00482                IF( NOUNIT ) THEN
00483                   D11 = A( KC )
00484                   D22 = A( KCNEXT )
00485                   D21 = A( KC+1 )
00486                   D12 = D21
00487                   DO 110 J = 1, NRHS
00488                      T1 = B( K, J )
00489                      T2 = B( K+1, J )
00490                      B( K, J ) = D11*T1 + D12*T2
00491                      B( K+1, J ) = D21*T1 + D22*T2
00492   110             CONTINUE
00493                END IF
00494                KC = KCNEXT + ( N-K )
00495                K = K + 2
00496             END IF
00497             GO TO 100
00498   120       CONTINUE
00499          END IF
00500 *
00501       END IF
00502       RETURN
00503 *
00504 *     End of ZLAVSP
00505 *
00506       END