File:  [local] / rpl / lapack / blas / dsyr2k.f
Revision 1.1.1.1 (vendor branch): download - view: text, annotated - select for diffs - revision graph
Tue Jan 26 15:22:45 2010 UTC (14 years, 9 months ago) by bertrand
Branches: JKB
CVS tags: start, rpl-4_0_14, rpl-4_0_13, rpl-4_0_12, rpl-4_0_11, rpl-4_0_10


Commit initial.

    1:       SUBROUTINE DSYR2K(UPLO,TRANS,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
    2: *     .. Scalar Arguments ..
    3:       DOUBLE PRECISION ALPHA,BETA
    4:       INTEGER K,LDA,LDB,LDC,N
    5:       CHARACTER TRANS,UPLO
    6: *     ..
    7: *     .. Array Arguments ..
    8:       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
    9: *     ..
   10: *
   11: *  Purpose
   12: *  =======
   13: *
   14: *  DSYR2K  performs one of the symmetric rank 2k operations
   15: *
   16: *     C := alpha*A*B' + alpha*B*A' + beta*C,
   17: *
   18: *  or
   19: *
   20: *     C := alpha*A'*B + alpha*B'*A + beta*C,
   21: *
   22: *  where  alpha and beta  are scalars, C is an  n by n  symmetric matrix
   23: *  and  A and B  are  n by k  matrices  in the  first  case  and  k by n
   24: *  matrices in the second case.
   25: *
   26: *  Arguments
   27: *  ==========
   28: *
   29: *  UPLO   - CHARACTER*1.
   30: *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
   31: *           triangular  part  of the  array  C  is to be  referenced  as
   32: *           follows:
   33: *
   34: *              UPLO = 'U' or 'u'   Only the  upper triangular part of  C
   35: *                                  is to be referenced.
   36: *
   37: *              UPLO = 'L' or 'l'   Only the  lower triangular part of  C
   38: *                                  is to be referenced.
   39: *
   40: *           Unchanged on exit.
   41: *
   42: *  TRANS  - CHARACTER*1.
   43: *           On entry,  TRANS  specifies the operation to be performed as
   44: *           follows:
   45: *
   46: *              TRANS = 'N' or 'n'   C := alpha*A*B' + alpha*B*A' +
   47: *                                        beta*C.
   48: *
   49: *              TRANS = 'T' or 't'   C := alpha*A'*B + alpha*B'*A +
   50: *                                        beta*C.
   51: *
   52: *              TRANS = 'C' or 'c'   C := alpha*A'*B + alpha*B'*A +
   53: *                                        beta*C.
   54: *
   55: *           Unchanged on exit.
   56: *
   57: *  N      - INTEGER.
   58: *           On entry,  N specifies the order of the matrix C.  N must be
   59: *           at least zero.
   60: *           Unchanged on exit.
   61: *
   62: *  K      - INTEGER.
   63: *           On entry with  TRANS = 'N' or 'n',  K  specifies  the number
   64: *           of  columns  of the  matrices  A and B,  and on  entry  with
   65: *           TRANS = 'T' or 't' or 'C' or 'c',  K  specifies  the  number
   66: *           of rows of the matrices  A and B.  K must be at least  zero.
   67: *           Unchanged on exit.
   68: *
   69: *  ALPHA  - DOUBLE PRECISION.
   70: *           On entry, ALPHA specifies the scalar alpha.
   71: *           Unchanged on exit.
   72: *
   73: *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
   74: *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
   75: *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
   76: *           part of the array  A  must contain the matrix  A,  otherwise
   77: *           the leading  k by n  part of the array  A  must contain  the
   78: *           matrix A.
   79: *           Unchanged on exit.
   80: *
   81: *  LDA    - INTEGER.
   82: *           On entry, LDA specifies the first dimension of A as declared
   83: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
   84: *           then  LDA must be at least  max( 1, n ), otherwise  LDA must
   85: *           be at least  max( 1, k ).
   86: *           Unchanged on exit.
   87: *
   88: *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
   89: *           k  when  TRANS = 'N' or 'n',  and is  n  otherwise.
   90: *           Before entry with  TRANS = 'N' or 'n',  the  leading  n by k
   91: *           part of the array  B  must contain the matrix  B,  otherwise
   92: *           the leading  k by n  part of the array  B  must contain  the
   93: *           matrix B.
   94: *           Unchanged on exit.
   95: *
   96: *  LDB    - INTEGER.
   97: *           On entry, LDB specifies the first dimension of B as declared
   98: *           in  the  calling  (sub)  program.   When  TRANS = 'N' or 'n'
   99: *           then  LDB must be at least  max( 1, n ), otherwise  LDB must
  100: *           be at least  max( 1, k ).
  101: *           Unchanged on exit.
  102: *
  103: *  BETA   - DOUBLE PRECISION.
  104: *           On entry, BETA specifies the scalar beta.
  105: *           Unchanged on exit.
  106: *
  107: *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
  108: *           Before entry  with  UPLO = 'U' or 'u',  the leading  n by n
  109: *           upper triangular part of the array C must contain the upper
  110: *           triangular part  of the  symmetric matrix  and the strictly
  111: *           lower triangular part of C is not referenced.  On exit, the
  112: *           upper triangular part of the array  C is overwritten by the
  113: *           upper triangular part of the updated matrix.
  114: *           Before entry  with  UPLO = 'L' or 'l',  the leading  n by n
  115: *           lower triangular part of the array C must contain the lower
  116: *           triangular part  of the  symmetric matrix  and the strictly
  117: *           upper triangular part of C is not referenced.  On exit, the
  118: *           lower triangular part of the array  C is overwritten by the
  119: *           lower triangular part of the updated matrix.
  120: *
  121: *  LDC    - INTEGER.
  122: *           On entry, LDC specifies the first dimension of C as declared
  123: *           in  the  calling  (sub)  program.   LDC  must  be  at  least
  124: *           max( 1, n ).
  125: *           Unchanged on exit.
  126: *
  127: *  Further Details
  128: *  ===============
  129: *
  130: *  Level 3 Blas routine.
  131: *
  132: *
  133: *  -- Written on 8-February-1989.
  134: *     Jack Dongarra, Argonne National Laboratory.
  135: *     Iain Duff, AERE Harwell.
  136: *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  137: *     Sven Hammarling, Numerical Algorithms Group Ltd.
  138: *
  139: *  =====================================================================
  140: *
  141: *     .. External Functions ..
  142:       LOGICAL LSAME
  143:       EXTERNAL LSAME
  144: *     ..
  145: *     .. External Subroutines ..
  146:       EXTERNAL XERBLA
  147: *     ..
  148: *     .. Intrinsic Functions ..
  149:       INTRINSIC MAX
  150: *     ..
  151: *     .. Local Scalars ..
  152:       DOUBLE PRECISION TEMP1,TEMP2
  153:       INTEGER I,INFO,J,L,NROWA
  154:       LOGICAL UPPER
  155: *     ..
  156: *     .. Parameters ..
  157:       DOUBLE PRECISION ONE,ZERO
  158:       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
  159: *     ..
  160: *
  161: *     Test the input parameters.
  162: *
  163:       IF (LSAME(TRANS,'N')) THEN
  164:           NROWA = N
  165:       ELSE
  166:           NROWA = K
  167:       END IF
  168:       UPPER = LSAME(UPLO,'U')
  169: *
  170:       INFO = 0
  171:       IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
  172:           INFO = 1
  173:       ELSE IF ((.NOT.LSAME(TRANS,'N')) .AND.
  174:      +         (.NOT.LSAME(TRANS,'T')) .AND.
  175:      +         (.NOT.LSAME(TRANS,'C'))) THEN
  176:           INFO = 2
  177:       ELSE IF (N.LT.0) THEN
  178:           INFO = 3
  179:       ELSE IF (K.LT.0) THEN
  180:           INFO = 4
  181:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
  182:           INFO = 7
  183:       ELSE IF (LDB.LT.MAX(1,NROWA)) THEN
  184:           INFO = 9
  185:       ELSE IF (LDC.LT.MAX(1,N)) THEN
  186:           INFO = 12
  187:       END IF
  188:       IF (INFO.NE.0) THEN
  189:           CALL XERBLA('DSYR2K',INFO)
  190:           RETURN
  191:       END IF
  192: *
  193: *     Quick return if possible.
  194: *
  195:       IF ((N.EQ.0) .OR. (((ALPHA.EQ.ZERO).OR.
  196:      +    (K.EQ.0)).AND. (BETA.EQ.ONE))) RETURN
  197: *
  198: *     And when  alpha.eq.zero.
  199: *
  200:       IF (ALPHA.EQ.ZERO) THEN
  201:           IF (UPPER) THEN
  202:               IF (BETA.EQ.ZERO) THEN
  203:                   DO 20 J = 1,N
  204:                       DO 10 I = 1,J
  205:                           C(I,J) = ZERO
  206:    10                 CONTINUE
  207:    20             CONTINUE
  208:               ELSE
  209:                   DO 40 J = 1,N
  210:                       DO 30 I = 1,J
  211:                           C(I,J) = BETA*C(I,J)
  212:    30                 CONTINUE
  213:    40             CONTINUE
  214:               END IF
  215:           ELSE
  216:               IF (BETA.EQ.ZERO) THEN
  217:                   DO 60 J = 1,N
  218:                       DO 50 I = J,N
  219:                           C(I,J) = ZERO
  220:    50                 CONTINUE
  221:    60             CONTINUE
  222:               ELSE
  223:                   DO 80 J = 1,N
  224:                       DO 70 I = J,N
  225:                           C(I,J) = BETA*C(I,J)
  226:    70                 CONTINUE
  227:    80             CONTINUE
  228:               END IF
  229:           END IF
  230:           RETURN
  231:       END IF
  232: *
  233: *     Start the operations.
  234: *
  235:       IF (LSAME(TRANS,'N')) THEN
  236: *
  237: *        Form  C := alpha*A*B' + alpha*B*A' + C.
  238: *
  239:           IF (UPPER) THEN
  240:               DO 130 J = 1,N
  241:                   IF (BETA.EQ.ZERO) THEN
  242:                       DO 90 I = 1,J
  243:                           C(I,J) = ZERO
  244:    90                 CONTINUE
  245:                   ELSE IF (BETA.NE.ONE) THEN
  246:                       DO 100 I = 1,J
  247:                           C(I,J) = BETA*C(I,J)
  248:   100                 CONTINUE
  249:                   END IF
  250:                   DO 120 L = 1,K
  251:                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
  252:                           TEMP1 = ALPHA*B(J,L)
  253:                           TEMP2 = ALPHA*A(J,L)
  254:                           DO 110 I = 1,J
  255:                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
  256:      +                                 B(I,L)*TEMP2
  257:   110                     CONTINUE
  258:                       END IF
  259:   120             CONTINUE
  260:   130         CONTINUE
  261:           ELSE
  262:               DO 180 J = 1,N
  263:                   IF (BETA.EQ.ZERO) THEN
  264:                       DO 140 I = J,N
  265:                           C(I,J) = ZERO
  266:   140                 CONTINUE
  267:                   ELSE IF (BETA.NE.ONE) THEN
  268:                       DO 150 I = J,N
  269:                           C(I,J) = BETA*C(I,J)
  270:   150                 CONTINUE
  271:                   END IF
  272:                   DO 170 L = 1,K
  273:                       IF ((A(J,L).NE.ZERO) .OR. (B(J,L).NE.ZERO)) THEN
  274:                           TEMP1 = ALPHA*B(J,L)
  275:                           TEMP2 = ALPHA*A(J,L)
  276:                           DO 160 I = J,N
  277:                               C(I,J) = C(I,J) + A(I,L)*TEMP1 +
  278:      +                                 B(I,L)*TEMP2
  279:   160                     CONTINUE
  280:                       END IF
  281:   170             CONTINUE
  282:   180         CONTINUE
  283:           END IF
  284:       ELSE
  285: *
  286: *        Form  C := alpha*A'*B + alpha*B'*A + C.
  287: *
  288:           IF (UPPER) THEN
  289:               DO 210 J = 1,N
  290:                   DO 200 I = 1,J
  291:                       TEMP1 = ZERO
  292:                       TEMP2 = ZERO
  293:                       DO 190 L = 1,K
  294:                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
  295:                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
  296:   190                 CONTINUE
  297:                       IF (BETA.EQ.ZERO) THEN
  298:                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
  299:                       ELSE
  300:                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
  301:      +                             ALPHA*TEMP2
  302:                       END IF
  303:   200             CONTINUE
  304:   210         CONTINUE
  305:           ELSE
  306:               DO 240 J = 1,N
  307:                   DO 230 I = J,N
  308:                       TEMP1 = ZERO
  309:                       TEMP2 = ZERO
  310:                       DO 220 L = 1,K
  311:                           TEMP1 = TEMP1 + A(L,I)*B(L,J)
  312:                           TEMP2 = TEMP2 + B(L,I)*A(L,J)
  313:   220                 CONTINUE
  314:                       IF (BETA.EQ.ZERO) THEN
  315:                           C(I,J) = ALPHA*TEMP1 + ALPHA*TEMP2
  316:                       ELSE
  317:                           C(I,J) = BETA*C(I,J) + ALPHA*TEMP1 +
  318:      +                             ALPHA*TEMP2
  319:                       END IF
  320:   230             CONTINUE
  321:   240         CONTINUE
  322:           END IF
  323:       END IF
  324: *
  325:       RETURN
  326: *
  327: *     End of DSYR2K.
  328: *
  329:       END

CVSweb interface <joel.bertrand@systella.fr>