Annotation of rpl/lapack/blas/zherk.f, revision 1.1.1.1

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

CVSweb interface <joel.bertrand@systella.fr>