Return to dsyconv.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1:       SUBROUTINE DSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
    2: *
    3: *  -- LAPACK PROTOTYPE routine (version 3.3.0) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2010
    7: *
    8: *  -- Written by Julie Langou of the Univ. of TN    --
    9: *
   10: *     .. Scalar Arguments ..
   11:       CHARACTER          UPLO, WAY
   12:       INTEGER            INFO, LDA, N
   13: *     ..
   14: *     .. Array Arguments ..
   15:       INTEGER            IPIV( * )
   16:       DOUBLE PRECISION   A( LDA, * ), WORK( * )
   17: *     ..
   18: *
   19: *  Purpose
   20: *  =======
   21: *
   22: *  DSYCONV convert A given by TRF into L and D and vice-versa.
   23: *  Get Non-diag elements of D (returned in workspace) and 
   24: *  apply or reverse permutation done in TRF.
   25: *
   26: *  Arguments
   27: *  =========
   28: *
   29: *  UPLO    (input) CHARACTER*1
   30: *          Specifies whether the details of the factorization are stored
   31: *          as an upper or lower triangular matrix.
   32: *          = 'U':  Upper triangular, form is A = U*D*U**T;
   33: *          = 'L':  Lower triangular, form is A = L*D*L**T.
   34: * 
   35: *  WAY     (input) CHARACTER*1
   36: *          = 'C': Convert 
   37: *          = 'R': Revert
   38: *
   39: *  N       (input) INTEGER
   40: *          The order of the matrix A.  N >= 0.
   41: *
   42: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
   43: *          The block diagonal matrix D and the multipliers used to
   44: *          obtain the factor U or L as computed by DSYTRF.
   45: *
   46: *  LDA     (input) INTEGER
   47: *          The leading dimension of the array A.  LDA >= max(1,N).
   48: *
   49: *  IPIV    (input) INTEGER array, dimension (N)
   50: *          Details of the interchanges and the block structure of D
   51: *          as determined by DSYTRF.
   52: *
   53: * WORK     (workspace) DOUBLE PRECISION array, dimension (N)
   54: *
   55: * LWORK    (input) INTEGER
   56: *          The length of WORK.  LWORK >=1. 
   57: *          LWORK = N
   58: *
   59: *          If LWORK = -1, then a workspace query is assumed; the routine
   60: *          only calculates the optimal size of the WORK array, returns
   61: *          this value as the first entry of the WORK array, and no error
   62: *          message related to LWORK is issued by XERBLA.
   63: *
   64: *  INFO    (output) INTEGER
   65: *          = 0:  successful exit
   66: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   67: *
   68: *  =====================================================================
   69: *
   70: *     .. Parameters ..
   71:       DOUBLE PRECISION   ZERO
   72:       PARAMETER          ( ZERO = 0.0D+0 )
   73: *     ..
   74: *     .. External Functions ..
   75:       LOGICAL            LSAME
   76:       EXTERNAL           LSAME
   77: *
   78: *     .. External Subroutines ..
   79:       EXTERNAL           XERBLA
   80: *     .. Local Scalars ..
   81:       LOGICAL            UPPER, CONVERT
   82:       INTEGER            I, IP, J
   83:       DOUBLE PRECISION   TEMP
   84: *     ..
   85: *     .. Executable Statements ..
   86: *
   87:       INFO = 0
   88:       UPPER = LSAME( UPLO, 'U' )
   89:       CONVERT = LSAME( WAY, 'C' )
   90:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
   91:          INFO = -1
   92:       ELSE IF( .NOT.CONVERT .AND. .NOT.LSAME( WAY, 'R' ) ) THEN
   93:          INFO = -2
   94:       ELSE IF( N.LT.0 ) THEN
   95:          INFO = -3
   96:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
   97:          INFO = -5
   98: 
   99:       END IF
  100:       IF( INFO.NE.0 ) THEN
  101:          CALL XERBLA( 'DSYCONV', -INFO )
  102:          RETURN
  103:       END IF
  104: *
  105: *     Quick return if possible
  106: *
  107:       IF( N.EQ.0 )
  108:      $   RETURN
  109: *
  110:       IF( UPPER ) THEN
  111: *
  112: *      A is UPPER
  113: *
  114: *      Convert A (A is upper)
  115: *
  116: *        Convert VALUE
  117: *
  118:          IF ( CONVERT ) THEN
  119:             I=N
  120:             WORK(1)=ZERO
  121:             DO WHILE ( I .GT. 1 )
  122:                IF( IPIV(I) .LT. 0 ) THEN
  123:                   WORK(I)=A(I-1,I)
  124:                   A(I-1,I)=ZERO
  125:                   I=I-1
  126:                ELSE
  127:                   WORK(I)=ZERO
  128:                ENDIF
  129:                I=I-1
  130:             END DO
  131: *
  132: *        Convert PERMUTATIONS
  133: *  
  134:          I=N
  135:          DO WHILE ( I .GE. 1 )
  136:             IF( IPIV(I) .GT. 0) THEN
  137:                IP=IPIV(I)
  138:                IF( I .LT. N) THEN
  139:                   DO 12 J= I+1,N
  140:                     TEMP=A(IP,J)
  141:                     A(IP,J)=A(I,J)
  142:                     A(I,J)=TEMP
  143:  12            CONTINUE
  144:                ENDIF
  145:             ELSE
  146:               IP=-IPIV(I)
  147:                IF( I .LT. N) THEN
  148:              DO 13 J= I+1,N
  149:                  TEMP=A(IP,J)
  150:                  A(IP,J)=A(I-1,J)
  151:                  A(I-1,J)=TEMP
  152:  13            CONTINUE
  153:                 ENDIF
  154:                 I=I-1
  155:            ENDIF
  156:            I=I-1
  157:         END DO
  158: 
  159:          ELSE
  160: *
  161: *      Revert A (A is upper)
  162: *
  163: *
  164: *        Revert PERMUTATIONS
  165: *  
  166:             I=1
  167:             DO WHILE ( I .LE. N )
  168:                IF( IPIV(I) .GT. 0 ) THEN
  169:                   IP=IPIV(I)
  170:                   IF( I .LT. N) THEN
  171:                   DO J= I+1,N
  172:                     TEMP=A(IP,J)
  173:                     A(IP,J)=A(I,J)
  174:                     A(I,J)=TEMP
  175:                   END DO
  176:                   ENDIF
  177:                ELSE
  178:                  IP=-IPIV(I)
  179:                  I=I+1
  180:                  IF( I .LT. N) THEN
  181:                     DO J= I+1,N
  182:                        TEMP=A(IP,J)
  183:                        A(IP,J)=A(I-1,J)
  184:                        A(I-1,J)=TEMP
  185:                     END DO
  186:                  ENDIF
  187:                ENDIF
  188:                I=I+1
  189:             END DO
  190: *
  191: *        Revert VALUE
  192: *
  193:             I=N
  194:             DO WHILE ( I .GT. 1 )
  195:                IF( IPIV(I) .LT. 0 ) THEN
  196:                   A(I-1,I)=WORK(I)
  197:                   I=I-1
  198:                ENDIF
  199:                I=I-1
  200:             END DO
  201:          END IF
  202:       ELSE
  203: *
  204: *      A is LOWER
  205: *
  206:          IF ( CONVERT ) THEN
  207: *
  208: *      Convert A (A is lower)
  209: *
  210: *
  211: *        Convert VALUE
  212: *
  213:             I=1
  214:             WORK(N)=ZERO
  215:             DO WHILE ( I .LE. N )
  216:                IF( I.LT.N .AND. IPIV(I) .LT. 0 ) THEN
  217:                   WORK(I)=A(I+1,I)
  218:                   A(I+1,I)=ZERO
  219:                   I=I+1
  220:                ELSE
  221:                   WORK(I)=ZERO
  222:                ENDIF
  223:                I=I+1
  224:             END DO
  225: *
  226: *        Convert PERMUTATIONS
  227: *
  228:          I=1
  229:          DO WHILE ( I .LE. N )
  230:             IF( IPIV(I) .GT. 0 ) THEN
  231:                IP=IPIV(I)
  232:                IF (I .GT. 1) THEN
  233:                DO 22 J= 1,I-1
  234:                  TEMP=A(IP,J)
  235:                  A(IP,J)=A(I,J)
  236:                  A(I,J)=TEMP
  237:  22            CONTINUE
  238:                ENDIF
  239:             ELSE
  240:               IP=-IPIV(I)
  241:               IF (I .GT. 1) THEN
  242:               DO 23 J= 1,I-1
  243:                  TEMP=A(IP,J)
  244:                  A(IP,J)=A(I+1,J)
  245:                  A(I+1,J)=TEMP
  246:  23           CONTINUE
  247:               ENDIF
  248:               I=I+1
  249:            ENDIF
  250:            I=I+1
  251:         END DO
  252:          ELSE
  253: *
  254: *      Revert A (A is lower)
  255: *
  256: *
  257: *        Revert PERMUTATIONS
  258: *
  259:             I=N
  260:             DO WHILE ( I .GE. 1 )
  261:                IF( IPIV(I) .GT. 0 ) THEN
  262:                   IP=IPIV(I)
  263:                   IF (I .GT. 1) THEN
  264:                      DO J= 1,I-1
  265:                         TEMP=A(I,J)
  266:                         A(I,J)=A(IP,J)
  267:                         A(IP,J)=TEMP
  268:                      END DO
  269:                   ENDIF
  270:                ELSE
  271:                   IP=-IPIV(I)
  272:                   I=I-1
  273:                   IF (I .GT. 1) THEN
  274:                      DO J= 1,I-1
  275:                         TEMP=A(I+1,J)
  276:                         A(I+1,J)=A(IP,J)
  277:                         A(IP,J)=TEMP
  278:                      END DO
  279:                   ENDIF
  280:                ENDIF
  281:                I=I-1
  282:             END DO
  283: *
  284: *        Revert VALUE
  285: *
  286:             I=1
  287:             DO WHILE ( I .LE. N-1 )
  288:                IF( IPIV(I) .LT. ZERO ) THEN
  289:                   A(I+1,I)=WORK(I)
  290:                   I=I+1
  291:                ENDIF
  292:                I=I+1
  293:             END DO
  294:          END IF
  295:       END IF
  296: 
  297:       RETURN
  298: *
  299: *     End of DSYCONV
  300: *
  301:       END