The R Project SVN R

Rev

Rev 23175 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 23175 Rev 23236
Line 7121... Line 7121...
7121
*     RETURN
7121
*     RETURN
7122
*
7122
*
7123
*     End of LSAME
7123
*     End of LSAME
7124
*
7124
*
7125
      END
7125
      END
7126
      SUBROUTINE ZGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB,
-
 
7127
     $                   BETA, C, LDC )
-
 
7128
*     .. Scalar Arguments ..
-
 
7129
      CHARACTER*1        TRANSA, TRANSB
-
 
7130
      INTEGER            M, N, K, LDA, LDB, LDC
-
 
7131
      COMPLEX*16         ALPHA, BETA
-
 
7132
*     .. Array Arguments ..
-
 
7133
      COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
-
 
7134
*     ..
-
 
7135
*
-
 
7136
*  Purpose
-
 
7137
*  =======
-
 
7138
*
-
 
7139
*  ZGEMM  performs one of the matrix-matrix operations
-
 
7140
*
-
 
7141
*     C := alpha*op( A )*op( B ) + beta*C,
-
 
7142
*
-
 
7143
*  where  op( X ) is one of
-
 
7144
*
-
 
7145
*     op( X ) = X   or   op( X ) = X'   or   op( X ) = conjg( X' ),
-
 
7146
*
-
 
7147
*  alpha and beta are scalars, and A, B and C are matrices, with op( A )
-
 
7148
*  an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.
-
 
7149
*
-
 
7150
*  Parameters
-
 
7151
*  ==========
-
 
7152
*
-
 
7153
*  TRANSA - CHARACTER*1.
-
 
7154
*           On entry, TRANSA specifies the form of op( A ) to be used in
-
 
7155
*           the matrix multiplication as follows:
-
 
7156
*
-
 
7157
*              TRANSA = 'N' or 'n',  op( A ) = A.
-
 
7158
*
-
 
7159
*              TRANSA = 'T' or 't',  op( A ) = A'.
-
 
7160
*
-
 
7161
*              TRANSA = 'C' or 'c',  op( A ) = conjg( A' ).
-
 
7162
*
-
 
7163
*           Unchanged on exit.
-
 
7164
*
-
 
7165
*  TRANSB - CHARACTER*1.
-
 
7166
*           On entry, TRANSB specifies the form of op( B ) to be used in
-
 
7167
*           the matrix multiplication as follows:
-
 
7168
*
-
 
7169
*              TRANSB = 'N' or 'n',  op( B ) = B.
-
 
7170
*
-
 
7171
*              TRANSB = 'T' or 't',  op( B ) = B'.
-
 
7172
*
-
 
7173
*              TRANSB = 'C' or 'c',  op( B ) = conjg( B' ).
-
 
7174
*
-
 
7175
*           Unchanged on exit.
-
 
7176
*
-
 
7177
*  M      - INTEGER.
-
 
7178
*           On entry,  M  specifies  the number  of rows  of the  matrix
-
 
7179
*           op( A )  and of the  matrix  C.  M  must  be at least  zero.
-
 
7180
*           Unchanged on exit.
-
 
7181
*
-
 
7182
*  N      - INTEGER.
-
 
7183
*           On entry,  N  specifies the number  of columns of the matrix
-
 
7184
*           op( B ) and the number of columns of the matrix C. N must be
-
 
7185
*           at least zero.
-
 
7186
*           Unchanged on exit.
-
 
7187
*
-
 
7188
*  K      - INTEGER.
-
 
7189
*           On entry,  K  specifies  the number of columns of the matrix
-
 
7190
*           op( A ) and the number of rows of the matrix op( B ). K must
-
 
7191
*           be at least  zero.
-
 
7192
*           Unchanged on exit.
-
 
7193
*
-
 
7194
*  ALPHA  - COMPLEX*16      .
-
 
7195
*           On entry, ALPHA specifies the scalar alpha.
-
 
7196
*           Unchanged on exit.
-
 
7197
*
-
 
7198
*  A      - COMPLEX*16       array of DIMENSION ( LDA, ka ), where ka is
-
 
7199
*           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
-
 
7200
*           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
-
 
7201
*           part of the array  A  must contain the matrix  A,  otherwise
-
 
7202
*           the leading  k by m  part of the array  A  must contain  the
-
 
7203
*           matrix A.
-
 
7204
*           Unchanged on exit.
-
 
7205
*
-
 
7206
*  LDA    - INTEGER.
-
 
7207
*           On entry, LDA specifies the first dimension of A as declared
-
 
7208
*           in the calling (sub) program. When  TRANSA = 'N' or 'n' then
-
 
7209
*           LDA must be at least  max( 1, m ), otherwise  LDA must be at
-
 
7210
*           least  max( 1, k ).
-
 
7211
*           Unchanged on exit.
-
 
7212
*
-
 
7213
*  B      - COMPLEX*16       array of DIMENSION ( LDB, kb ), where kb is
-
 
7214
*           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
-
 
7215
*           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
-
 
7216
*           part of the array  B  must contain the matrix  B,  otherwise
-
 
7217
*           the leading  n by k  part of the array  B  must contain  the
-
 
7218
*           matrix B.
-
 
7219
*           Unchanged on exit.
-
 
7220
*
-
 
7221
*  LDB    - INTEGER.
-
 
7222
*           On entry, LDB specifies the first dimension of B as declared
-
 
7223
*           in the calling (sub) program. When  TRANSB = 'N' or 'n' then
-
 
7224
*           LDB must be at least  max( 1, k ), otherwise  LDB must be at
-
 
7225
*           least  max( 1, n ).
-
 
7226
*           Unchanged on exit.
-
 
7227
*
-
 
7228
*  BETA   - COMPLEX*16      .
-
 
7229
*           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
-
 
7230
*           supplied as zero then C need not be set on input.
-
 
7231
*           Unchanged on exit.
-
 
7232
*
-
 
7233
*  C      - COMPLEX*16       array of DIMENSION ( LDC, n ).
-
 
7234
*           Before entry, the leading  m by n  part of the array  C must
-
 
7235
*           contain the matrix  C,  except when  beta  is zero, in which
-
 
7236
*           case C need not be set on entry.
-
 
7237
*           On exit, the array  C  is overwritten by the  m by n  matrix
-
 
7238
*           ( alpha*op( A )*op( B ) + beta*C ).
-
 
7239
*
-
 
7240
*  LDC    - INTEGER.
-
 
7241
*           On entry, LDC specifies the first dimension of C as declared
-
 
7242
*           in  the  calling  (sub)  program.   LDC  must  be  at  least
-
 
7243
*           max( 1, m ).
-
 
7244
*           Unchanged on exit.
-
 
7245
*
-
 
7246
*
-
 
7247
*  Level 3 Blas routine.
-
 
7248
*
-
 
7249
*  -- Written on 8-February-1989.
-
 
7250
*     Jack Dongarra, Argonne National Laboratory.
-
 
7251
*     Iain Duff, AERE Harwell.
-
 
7252
*     Jeremy Du Croz, Numerical Algorithms Group Ltd.
-
 
7253
*     Sven Hammarling, Numerical Algorithms Group Ltd.
-
 
7254
*
-
 
7255
*
-
 
7256
*     .. External Functions ..
-
 
7257
      LOGICAL            LSAME
-
 
7258
      EXTERNAL           LSAME
-
 
7259
*     .. External Subroutines ..
-
 
7260
      EXTERNAL           XERBLA
-
 
7261
*     .. Intrinsic Functions ..
-
 
7262
      INTRINSIC          DCONJG, MAX
-
 
7263
*     .. Local Scalars ..
-
 
7264
      LOGICAL            CONJA, CONJB, NOTA, NOTB
-
 
7265
      INTEGER            I, INFO, J, L, NCOLA, NROWA, NROWB
-
 
7266
      COMPLEX*16         TEMP
-
 
7267
*     .. Parameters ..
-
 
7268
      COMPLEX*16         ONE
-
 
7269
      PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
-
 
7270
      COMPLEX*16         ZERO
-
 
7271
      PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
-
 
7272
*     ..
-
 
7273
*     .. Executable Statements ..
-
 
7274
*
-
 
7275
*     Set  NOTA  and  NOTB  as  true if  A  and  B  respectively are not
-
 
7276
*     conjugated or transposed, set  CONJA and CONJB  as true if  A  and
-
 
7277
*     B  respectively are to be  transposed but  not conjugated  and set
-
 
7278
*     NROWA, NCOLA and  NROWB  as the number of rows and  columns  of  A
-
 
7279
*     and the number of rows of  B  respectively.
-
 
7280
*
-
 
7281
      NOTA  = LSAME( TRANSA, 'N' )
-
 
7282
      NOTB  = LSAME( TRANSB, 'N' )
-
 
7283
      CONJA = LSAME( TRANSA, 'C' )
-
 
7284
      CONJB = LSAME( TRANSB, 'C' )
-
 
7285
      IF( NOTA )THEN
-
 
7286
         NROWA = M
-
 
7287
         NCOLA = K
-
 
7288
      ELSE
-
 
7289
         NROWA = K
-
 
7290
         NCOLA = M
-
 
7291
      END IF
-
 
7292
      IF( NOTB )THEN
-
 
7293
         NROWB = K
-
 
7294
      ELSE
-
 
7295
         NROWB = N
-
 
7296
      END IF
-
 
7297
*
-
 
7298
*     Test the input parameters.
-
 
7299
*
-
 
7300
      INFO = 0
-
 
7301
      IF(      ( .NOT.NOTA                 ).AND.
-
 
7302
     $         ( .NOT.CONJA                ).AND.
-
 
7303
     $         ( .NOT.LSAME( TRANSA, 'T' ) )      )THEN
-
 
7304
         INFO = 1
-
 
7305
      ELSE IF( ( .NOT.NOTB                 ).AND.
-
 
7306
     $         ( .NOT.CONJB                ).AND.
-
 
7307
     $         ( .NOT.LSAME( TRANSB, 'T' ) )      )THEN
-
 
7308
         INFO = 2
-
 
7309
      ELSE IF( M  .LT.0               )THEN
-
 
7310
         INFO = 3
-
 
7311
      ELSE IF( N  .LT.0               )THEN
-
 
7312
         INFO = 4
-
 
7313
      ELSE IF( K  .LT.0               )THEN
-
 
7314
         INFO = 5
-
 
7315
      ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
-
 
7316
         INFO = 8
-
 
7317
      ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN
-
 
7318
         INFO = 10
-
 
7319
      ELSE IF( LDC.LT.MAX( 1, M     ) )THEN
-
 
7320
         INFO = 13
-
 
7321
      END IF
-
 
7322
      IF( INFO.NE.0 )THEN
-
 
7323
         CALL XERBLA( 'ZGEMM ', INFO )
-
 
7324
         RETURN
-
 
7325
      END IF
-
 
7326
*
-
 
7327
*     Quick return if possible.
-
 
7328
*
-
 
7329
      IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
-
 
7330
     $    ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) )
-
 
7331
     $   RETURN
-
 
7332
*
-
 
7333
*     And when  alpha.eq.zero.
-
 
7334
*
-
 
7335
      IF( ALPHA.EQ.ZERO )THEN
-
 
7336
         IF( BETA.EQ.ZERO )THEN
-
 
7337
            DO 20, J = 1, N
-
 
7338
               DO 10, I = 1, M
-
 
7339
                  C( I, J ) = ZERO
-
 
7340
   10          CONTINUE
-
 
7341
   20       CONTINUE
-
 
7342
         ELSE
-
 
7343
            DO 40, J = 1, N
-
 
7344
               DO 30, I = 1, M
-
 
7345
                  C( I, J ) = BETA*C( I, J )
-
 
7346
   30          CONTINUE
-
 
7347
   40       CONTINUE
-
 
7348
         END IF
-
 
7349
         RETURN
-
 
7350
      END IF
-
 
7351
*
-
 
7352
*     Start the operations.
-
 
7353
*
-
 
7354
      IF( NOTB )THEN
-
 
7355
         IF( NOTA )THEN
-
 
7356
*
-
 
7357
*           Form  C := alpha*A*B + beta*C.
-
 
7358
*
-
 
7359
            DO 90, J = 1, N
-
 
7360
               IF( BETA.EQ.ZERO )THEN
-
 
7361
                  DO 50, I = 1, M
-
 
7362
                     C( I, J ) = ZERO
-
 
7363
   50             CONTINUE
-
 
7364
               ELSE IF( BETA.NE.ONE )THEN
-
 
7365
                  DO 60, I = 1, M
-
 
7366
                     C( I, J ) = BETA*C( I, J )
-
 
7367
   60             CONTINUE
-
 
7368
               END IF
-
 
7369
               DO 80, L = 1, K
-
 
7370
                  IF( B( L, J ).NE.ZERO )THEN
-
 
7371
                     TEMP = ALPHA*B( L, J )
-
 
7372
                     DO 70, I = 1, M
-
 
7373
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
-
 
7374
   70                CONTINUE
-
 
7375
                  END IF
-
 
7376
   80          CONTINUE
-
 
7377
   90       CONTINUE
-
 
7378
         ELSE IF( CONJA )THEN
-
 
7379
*
-
 
7380
*           Form  C := alpha*conjg( A' )*B + beta*C.
-
 
7381
*
-
 
7382
            DO 120, J = 1, N
-
 
7383
               DO 110, I = 1, M
-
 
7384
                  TEMP = ZERO
-
 
7385
                  DO 100, L = 1, K
-
 
7386
                     TEMP = TEMP + DCONJG( A( L, I ) )*B( L, J )
-
 
7387
  100             CONTINUE
-
 
7388
                  IF( BETA.EQ.ZERO )THEN
-
 
7389
                     C( I, J ) = ALPHA*TEMP
-
 
7390
                  ELSE
-
 
7391
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7392
                  END IF
-
 
7393
  110          CONTINUE
-
 
7394
  120       CONTINUE
-
 
7395
         ELSE
-
 
7396
*
-
 
7397
*           Form  C := alpha*A'*B + beta*C
-
 
7398
*
-
 
7399
            DO 150, J = 1, N
-
 
7400
               DO 140, I = 1, M
-
 
7401
                  TEMP = ZERO
-
 
7402
                  DO 130, L = 1, K
-
 
7403
                     TEMP = TEMP + A( L, I )*B( L, J )
-
 
7404
  130             CONTINUE
-
 
7405
                  IF( BETA.EQ.ZERO )THEN
-
 
7406
                     C( I, J ) = ALPHA*TEMP
-
 
7407
                  ELSE
-
 
7408
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7409
                  END IF
-
 
7410
  140          CONTINUE
-
 
7411
  150       CONTINUE
-
 
7412
         END IF
-
 
7413
      ELSE IF( NOTA )THEN
-
 
7414
         IF( CONJB )THEN
-
 
7415
*
-
 
7416
*           Form  C := alpha*A*conjg( B' ) + beta*C.
-
 
7417
*
-
 
7418
            DO 200, J = 1, N
-
 
7419
               IF( BETA.EQ.ZERO )THEN
-
 
7420
                  DO 160, I = 1, M
-
 
7421
                     C( I, J ) = ZERO
-
 
7422
  160             CONTINUE
-
 
7423
               ELSE IF( BETA.NE.ONE )THEN
-
 
7424
                  DO 170, I = 1, M
-
 
7425
                     C( I, J ) = BETA*C( I, J )
-
 
7426
  170             CONTINUE
-
 
7427
               END IF
-
 
7428
               DO 190, L = 1, K
-
 
7429
                  IF( B( J, L ).NE.ZERO )THEN
-
 
7430
                     TEMP = ALPHA*DCONJG( B( J, L ) )
-
 
7431
                     DO 180, I = 1, M
-
 
7432
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
-
 
7433
  180                CONTINUE
-
 
7434
                  END IF
-
 
7435
  190          CONTINUE
-
 
7436
  200       CONTINUE
-
 
7437
         ELSE
-
 
7438
*
-
 
7439
*           Form  C := alpha*A*B'          + beta*C
-
 
7440
*
-
 
7441
            DO 250, J = 1, N
-
 
7442
               IF( BETA.EQ.ZERO )THEN
-
 
7443
                  DO 210, I = 1, M
-
 
7444
                     C( I, J ) = ZERO
-
 
7445
  210             CONTINUE
-
 
7446
               ELSE IF( BETA.NE.ONE )THEN
-
 
7447
                  DO 220, I = 1, M
-
 
7448
                     C( I, J ) = BETA*C( I, J )
-
 
7449
  220             CONTINUE
-
 
7450
               END IF
-
 
7451
               DO 240, L = 1, K
-
 
7452
                  IF( B( J, L ).NE.ZERO )THEN
-
 
7453
                     TEMP = ALPHA*B( J, L )
-
 
7454
                     DO 230, I = 1, M
-
 
7455
                        C( I, J ) = C( I, J ) + TEMP*A( I, L )
-
 
7456
  230                CONTINUE
-
 
7457
                  END IF
-
 
7458
  240          CONTINUE
-
 
7459
  250       CONTINUE
-
 
7460
         END IF
-
 
7461
      ELSE IF( CONJA )THEN
-
 
7462
         IF( CONJB )THEN
-
 
7463
*
-
 
7464
*           Form  C := alpha*conjg( A' )*conjg( B' ) + beta*C.
-
 
7465
*
-
 
7466
            DO 280, J = 1, N
-
 
7467
               DO 270, I = 1, M
-
 
7468
                  TEMP = ZERO
-
 
7469
                  DO 260, L = 1, K
-
 
7470
                     TEMP = TEMP +
-
 
7471
     $                      DCONJG( A( L, I ) )*DCONJG( B( J, L ) )
-
 
7472
  260             CONTINUE
-
 
7473
                  IF( BETA.EQ.ZERO )THEN
-
 
7474
                     C( I, J ) = ALPHA*TEMP
-
 
7475
                  ELSE
-
 
7476
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7477
                  END IF
-
 
7478
  270          CONTINUE
-
 
7479
  280       CONTINUE
-
 
7480
         ELSE
-
 
7481
*
-
 
7482
*           Form  C := alpha*conjg( A' )*B' + beta*C
-
 
7483
*
-
 
7484
            DO 310, J = 1, N
-
 
7485
               DO 300, I = 1, M
-
 
7486
                  TEMP = ZERO
-
 
7487
                  DO 290, L = 1, K
-
 
7488
                     TEMP = TEMP + DCONJG( A( L, I ) )*B( J, L )
-
 
7489
  290             CONTINUE
-
 
7490
                  IF( BETA.EQ.ZERO )THEN
-
 
7491
                     C( I, J ) = ALPHA*TEMP
-
 
7492
                  ELSE
-
 
7493
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7494
                  END IF
-
 
7495
  300          CONTINUE
-
 
7496
  310       CONTINUE
-
 
7497
         END IF
-
 
7498
      ELSE
-
 
7499
         IF( CONJB )THEN
-
 
7500
*
-
 
7501
*           Form  C := alpha*A'*conjg( B' ) + beta*C
-
 
7502
*
-
 
7503
            DO 340, J = 1, N
-
 
7504
               DO 330, I = 1, M
-
 
7505
                  TEMP = ZERO
-
 
7506
                  DO 320, L = 1, K
-
 
7507
                     TEMP = TEMP + A( L, I )*DCONJG( B( J, L ) )
-
 
7508
  320             CONTINUE
-
 
7509
                  IF( BETA.EQ.ZERO )THEN
-
 
7510
                     C( I, J ) = ALPHA*TEMP
-
 
7511
                  ELSE
-
 
7512
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7513
                  END IF
-
 
7514
  330          CONTINUE
-
 
7515
  340       CONTINUE
-
 
7516
         ELSE
-
 
7517
*
-
 
7518
*           Form  C := alpha*A'*B' + beta*C
-
 
7519
*
-
 
7520
            DO 370, J = 1, N
-
 
7521
               DO 360, I = 1, M
-
 
7522
                  TEMP = ZERO
-
 
7523
                  DO 350, L = 1, K
-
 
7524
                     TEMP = TEMP + A( L, I )*B( J, L )
-
 
7525
  350             CONTINUE
-
 
7526
                  IF( BETA.EQ.ZERO )THEN
-
 
7527
                     C( I, J ) = ALPHA*TEMP
-
 
7528
                  ELSE
-
 
7529
                     C( I, J ) = ALPHA*TEMP + BETA*C( I, J )
-
 
7530
                  END IF
-
 
7531
  360          CONTINUE
-
 
7532
  370       CONTINUE
-
 
7533
         END IF
-
 
7534
      END IF
-
 
7535
*
-
 
7536
      RETURN
-
 
7537
*
-
 
7538
*     End of ZGEMM .
-
 
7539
*
-
 
7540
      END
-