LAPACK: TESTING/LIN/clavsp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLAVSP( 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            A( * ), B( LDB, * )
00015 *     ..
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *     CLAVSP  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 CSPTRF.
00024 *     CSPTRF 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', CLAVSP multiplies either by U or U * D
00034 *     (or L or L * D).
00035 *     If TRANS = 'C' or 'c', CLAVSP 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 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 CSPTRF.
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 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            ONE
00110       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
00111 *     ..
00112 *     .. Local Scalars ..
00113       LOGICAL            NOUNIT
00114       INTEGER            J, K, KC, KCNEXT, KP
00115       COMPLEX            D11, D12, D21, D22, T1, T2
00116 *     ..
00117 *     .. External Functions ..
00118       LOGICAL            LSAME
00119       EXTERNAL           LSAME
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           CGEMV, CGERU, CSCAL, CSWAP, XERBLA
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( 'CLAVSP ', -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 CSCAL( 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 CGERU( 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 CSWAP( 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 CGERU( K-1, NRHS, ONE, A( KC ), 1, B( K, 1 ),
00230      $                        LDB, B( 1, 1 ), LDB )
00231                   CALL CGERU( 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 CSWAP( 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 CSCAL( 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 CGERU( 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 CSWAP( 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 CGERU( N-K, NRHS, ONE, A( KC+1 ), 1, B( K, 1 ),
00318      $                        LDB, B( K+1, 1 ), LDB )
00319                   CALL CGERU( 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 CSWAP( 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       IF( K.LT.1 )
00353      $         GO TO 90
00354             KC = KC - K
00355 *
00356 *           1 x 1 pivot block.
00357 *
00358             IF( IPIV( K ).GT.0 ) THEN
00359                IF( K.GT.1 ) THEN
00360 *
00361 *                 Interchange if P(K) != I.
00362 *
00363                   KP = IPIV( K )
00364                   IF( KP.NE.K )
00365      $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00366 *
00367 *                 Apply the transformation:
00368 *                    y := y - B' * conjg(x)
00369 *                 where x is a column of A and y is a row of B.
00370 *
00371                   CALL CGEMV( 'Transpose', K-1, NRHS, ONE, B, LDB,
00372      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00373                END IF
00374                IF( NOUNIT )
00375      $            CALL CSCAL( NRHS, A( KC+K-1 ), B( K, 1 ), LDB )
00376                K = K - 1
00377 *
00378 *           2 x 2 pivot block.
00379 *
00380             ELSE
00381                KCNEXT = KC - ( K-1 )
00382                IF( K.GT.2 ) THEN
00383 *
00384 *                 Interchange if P(K) != I.
00385 *
00386                   KP = ABS( IPIV( K ) )
00387                   IF( KP.NE.K-1 )
00388      $               CALL CSWAP( NRHS, B( K-1, 1 ), LDB, B( KP, 1 ),
00389      $                           LDB )
00390 *
00391 *                 Apply the transformations.
00392 *
00393                   CALL CGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00394      $                        A( KC ), 1, ONE, B( K, 1 ), LDB )
00395 *
00396                   CALL CGEMV( 'Transpose', K-2, NRHS, ONE, B, LDB,
00397      $                        A( KCNEXT ), 1, ONE, B( K-1, 1 ), LDB )
00398                END IF
00399 *
00400 *              Multiply by the diagonal block if non-unit.
00401 *
00402                IF( NOUNIT ) THEN
00403                   D11 = A( KC-1 )
00404                   D22 = A( KC+K-1 )
00405                   D12 = A( KC+K-2 )
00406                   D21 = D12
00407                   DO 80 J = 1, NRHS
00408                      T1 = B( K-1, J )
00409                      T2 = B( K, J )
00410                      B( K-1, J ) = D11*T1 + D12*T2
00411                      B( K, J ) = D21*T1 + D22*T2
00412    80             CONTINUE
00413                END IF
00414                KC = KCNEXT
00415                K = K - 2
00416             END IF
00417             GO TO 70
00418    90       CONTINUE
00419 *
00420 *        Form  B := L^T*B
00421 *        where L  = P(1)*inv(L(1))* ... *P(m)*inv(L(m))
00422 *        and   L^T = inv(L(m))*P(m)* ... *inv(L(1))*P(1)
00423 *
00424          ELSE
00425 *
00426 *           Loop forward applying the L-transformations.
00427 *
00428             K = 1
00429             KC = 1
00430   100       CONTINUE
00431             IF( K.GT.N )
00432      $         GO TO 120
00433 *
00434 *           1 x 1 pivot block
00435 *
00436             IF( IPIV( K ).GT.0 ) THEN
00437                IF( K.LT.N ) THEN
00438 *
00439 *                 Interchange if P(K) != I.
00440 *
00441                   KP = IPIV( K )
00442                   IF( KP.NE.K )
00443      $               CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
00444 *
00445 *                 Apply the transformation
00446 *
00447                   CALL CGEMV( 'Transpose', N-K, NRHS, ONE, B( K+1, 1 ),
00448      $                        LDB, A( KC+1 ), 1, ONE, B( K, 1 ), LDB )
00449                END IF
00450                IF( NOUNIT )
00451      $            CALL CSCAL( NRHS, A( KC ), B( K, 1 ), LDB )
00452                KC = KC + N - K + 1
00453                K = K + 1
00454 *
00455 *           2 x 2 pivot block.
00456 *
00457             ELSE
00458                KCNEXT = KC + N - K + 1
00459                IF( K.LT.N-1 ) THEN
00460 *
00461 *              Interchange if P(K) != I.
00462 *
00463                   KP = ABS( IPIV( K ) )
00464                   IF( KP.NE.K+1 )
00465      $               CALL CSWAP( NRHS, B( K+1, 1 ), LDB, B( KP, 1 ),
00466      $                           LDB )
00467 *
00468 *                 Apply the transformation
00469 *
00470                   CALL CGEMV( 'Transpose', N-K-1, NRHS, ONE,
00471      $                        B( K+2, 1 ), LDB, A( KCNEXT+1 ), 1, ONE,
00472      $                        B( K+1, 1 ), LDB )
00473 *
00474                   CALL CGEMV( 'Transpose', N-K-1, NRHS, ONE,
00475      $                        B( K+2, 1 ), LDB, A( KC+2 ), 1, ONE,
00476      $                        B( K, 1 ), LDB )
00477                END IF
00478 *
00479 *              Multiply by the diagonal block if non-unit.
00480 *
00481                IF( NOUNIT ) THEN
00482                   D11 = A( KC )
00483                   D22 = A( KCNEXT )
00484                   D21 = A( KC+1 )
00485                   D12 = D21
00486                   DO 110 J = 1, NRHS
00487                      T1 = B( K, J )
00488                      T2 = B( K+1, J )
00489                      B( K, J ) = D11*T1 + D12*T2
00490                      B( K+1, J ) = D21*T1 + D22*T2
00491   110             CONTINUE
00492                END IF
00493                KC = KCNEXT + ( N-K )
00494                K = K + 2
00495             END IF
00496             GO TO 100
00497   120       CONTINUE
00498          END IF
00499 *
00500       END IF
00501       RETURN
00502 *
00503 *     End of CLAVSP
00504 *
00505       END