LAPACK 3.3.1
Linear Algebra PACKage

dsymm.f

Go to the documentation of this file.
00001       SUBROUTINE DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
00002 *     .. Scalar Arguments ..
00003       DOUBLE PRECISION ALPHA,BETA
00004       INTEGER LDA,LDB,LDC,M,N
00005       CHARACTER SIDE,UPLO
00006 *     ..
00007 *     .. Array Arguments ..
00008       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
00009 *     ..
00010 *
00011 *  Purpose
00012 *  =======
00013 *
00014 *  DSYMM  performs one of the matrix-matrix operations
00015 *
00016 *     C := alpha*A*B + beta*C,
00017 *
00018 *  or
00019 *
00020 *     C := alpha*B*A + beta*C,
00021 *
00022 *  where alpha and beta are scalars,  A is a symmetric matrix and  B and
00023 *  C are  m by n matrices.
00024 *
00025 *  Arguments
00026 *  ==========
00027 *
00028 *  SIDE   - CHARACTER*1.
00029 *           On entry,  SIDE  specifies whether  the  symmetric matrix  A
00030 *           appears on the  left or right  in the  operation as follows:
00031 *
00032 *              SIDE = 'L' or 'l'   C := alpha*A*B + beta*C,
00033 *
00034 *              SIDE = 'R' or 'r'   C := alpha*B*A + beta*C,
00035 *
00036 *           Unchanged on exit.
00037 *
00038 *  UPLO   - CHARACTER*1.
00039 *           On  entry,   UPLO  specifies  whether  the  upper  or  lower
00040 *           triangular  part  of  the  symmetric  matrix   A  is  to  be
00041 *           referenced as follows:
00042 *
00043 *              UPLO = 'U' or 'u'   Only the upper triangular part of the
00044 *                                  symmetric matrix is to be referenced.
00045 *
00046 *              UPLO = 'L' or 'l'   Only the lower triangular part of the
00047 *                                  symmetric matrix is to be referenced.
00048 *
00049 *           Unchanged on exit.
00050 *
00051 *  M      - INTEGER.
00052 *           On entry,  M  specifies the number of rows of the matrix  C.
00053 *           M  must be at least zero.
00054 *           Unchanged on exit.
00055 *
00056 *  N      - INTEGER.
00057 *           On entry, N specifies the number of columns of the matrix C.
00058 *           N  must be at least zero.
00059 *           Unchanged on exit.
00060 *
00061 *  ALPHA  - DOUBLE PRECISION.
00062 *           On entry, ALPHA specifies the scalar alpha.
00063 *           Unchanged on exit.
00064 *
00065 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
00066 *           m  when  SIDE = 'L' or 'l'  and is  n otherwise.
00067 *           Before entry  with  SIDE = 'L' or 'l',  the  m by m  part of
00068 *           the array  A  must contain the  symmetric matrix,  such that
00069 *           when  UPLO = 'U' or 'u', the leading m by m upper triangular
00070 *           part of the array  A  must contain the upper triangular part
00071 *           of the  symmetric matrix and the  strictly  lower triangular
00072 *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
00073 *           the leading  m by m  lower triangular part  of the  array  A
00074 *           must  contain  the  lower triangular part  of the  symmetric
00075 *           matrix and the  strictly upper triangular part of  A  is not
00076 *           referenced.
00077 *           Before entry  with  SIDE = 'R' or 'r',  the  n by n  part of
00078 *           the array  A  must contain the  symmetric matrix,  such that
00079 *           when  UPLO = 'U' or 'u', the leading n by n upper triangular
00080 *           part of the array  A  must contain the upper triangular part
00081 *           of the  symmetric matrix and the  strictly  lower triangular
00082 *           part of  A  is not referenced,  and when  UPLO = 'L' or 'l',
00083 *           the leading  n by n  lower triangular part  of the  array  A
00084 *           must  contain  the  lower triangular part  of the  symmetric
00085 *           matrix and the  strictly upper triangular part of  A  is not
00086 *           referenced.
00087 *           Unchanged on exit.
00088 *
00089 *  LDA    - INTEGER.
00090 *           On entry, LDA specifies the first dimension of A as declared
00091 *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
00092 *           LDA must be at least  max( 1, m ), otherwise  LDA must be at
00093 *           least  max( 1, n ).
00094 *           Unchanged on exit.
00095 *
00096 *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
00097 *           Before entry, the leading  m by n part of the array  B  must
00098 *           contain the matrix B.
00099 *           Unchanged on exit.
00100 *
00101 *  LDB    - INTEGER.
00102 *           On entry, LDB specifies the first dimension of B as declared
00103 *           in  the  calling  (sub)  program.   LDB  must  be  at  least
00104 *           max( 1, m ).
00105 *           Unchanged on exit.
00106 *
00107 *  BETA   - DOUBLE PRECISION.
00108 *           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
00109 *           supplied as zero then C need not be set on input.
00110 *           Unchanged on exit.
00111 *
00112 *  C      - DOUBLE PRECISION array of DIMENSION ( LDC, n ).
00113 *           Before entry, the leading  m by n  part of the array  C must
00114 *           contain the matrix  C,  except when  beta  is zero, in which
00115 *           case C need not be set on entry.
00116 *           On exit, the array  C  is overwritten by the  m by n updated
00117 *           matrix.
00118 *
00119 *  LDC    - INTEGER.
00120 *           On entry, LDC specifies the first dimension of C as declared
00121 *           in  the  calling  (sub)  program.   LDC  must  be  at  least
00122 *           max( 1, m ).
00123 *           Unchanged on exit.
00124 *
00125 *  Further Details
00126 *  ===============
00127 *
00128 *  Level 3 Blas routine.
00129 *
00130 *  -- Written on 8-February-1989.
00131 *     Jack Dongarra, Argonne National Laboratory.
00132 *     Iain Duff, AERE Harwell.
00133 *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
00134 *     Sven Hammarling, Numerical Algorithms Group Ltd.
00135 *
00136 *  =====================================================================
00137 *
00138 *     .. External Functions ..
00139       LOGICAL LSAME
00140       EXTERNAL LSAME
00141 *     ..
00142 *     .. External Subroutines ..
00143       EXTERNAL XERBLA
00144 *     ..
00145 *     .. Intrinsic Functions ..
00146       INTRINSIC MAX
00147 *     ..
00148 *     .. Local Scalars ..
00149       DOUBLE PRECISION TEMP1,TEMP2
00150       INTEGER I,INFO,J,K,NROWA
00151       LOGICAL UPPER
00152 *     ..
00153 *     .. Parameters ..
00154       DOUBLE PRECISION ONE,ZERO
00155       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
00156 *     ..
00157 *
00158 *     Set NROWA as the number of rows of A.
00159 *
00160       IF (LSAME(SIDE,'L')) THEN
00161           NROWA = M
00162       ELSE
00163           NROWA = N
00164       END IF
00165       UPPER = LSAME(UPLO,'U')
00166 *
00167 *     Test the input parameters.
00168 *
00169       INFO = 0
00170       IF ((.NOT.LSAME(SIDE,'L')) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
00171           INFO = 1
00172       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
00173           INFO = 2
00174       ELSE IF (M.LT.0) THEN
00175           INFO = 3
00176       ELSE IF (N.LT.0) THEN
00177           INFO = 4
00178       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
00179           INFO = 7
00180       ELSE IF (LDB.LT.MAX(1,M)) THEN
00181           INFO = 9
00182       ELSE IF (LDC.LT.MAX(1,M)) THEN
00183           INFO = 12
00184       END IF
00185       IF (INFO.NE.0) THEN
00186           CALL XERBLA('DSYMM ',INFO)
00187           RETURN
00188       END IF
00189 *
00190 *     Quick return if possible.
00191 *
00192       IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
00193      +    ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
00194 *
00195 *     And when  alpha.eq.zero.
00196 *
00197       IF (ALPHA.EQ.ZERO) THEN
00198           IF (BETA.EQ.ZERO) THEN
00199               DO 20 J = 1,N
00200                   DO 10 I = 1,M
00201                       C(I,J) = ZERO
00202    10             CONTINUE
00203    20         CONTINUE
00204           ELSE
00205               DO 40 J = 1,N
00206                   DO 30 I = 1,M
00207                       C(I,J) = BETA*C(I,J)
00208    30             CONTINUE
00209    40         CONTINUE
00210           END IF
00211           RETURN
00212       END IF
00213 *
00214 *     Start the operations.
00215 *
00216       IF (LSAME(SIDE,'L')) THEN
00217 *
00218 *        Form  C := alpha*A*B + beta*C.
00219 *
00220           IF (UPPER) THEN
00221               DO 70 J = 1,N
00222                   DO 60 I = 1,M
00223                       TEMP1 = ALPHA*B(I,J)
00224                       TEMP2 = ZERO
00225                       DO 50 K = 1,I - 1
00226                           C(K,J) = C(K,J) + TEMP1*A(K,I)
00227                           TEMP2 = TEMP2 + B(K,J)*A(K,I)
00228    50                 CONTINUE
00229                       IF (BETA.EQ.ZERO) THEN
00230                           C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
00231                       ELSE
00232                           C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
00233      +                             ALPHA*TEMP2
00234                       END IF
00235    60             CONTINUE
00236    70         CONTINUE
00237           ELSE
00238               DO 100 J = 1,N
00239                   DO 90 I = M,1,-1
00240                       TEMP1 = ALPHA*B(I,J)
00241                       TEMP2 = ZERO
00242                       DO 80 K = I + 1,M
00243                           C(K,J) = C(K,J) + TEMP1*A(K,I)
00244                           TEMP2 = TEMP2 + B(K,J)*A(K,I)
00245    80                 CONTINUE
00246                       IF (BETA.EQ.ZERO) THEN
00247                           C(I,J) = TEMP1*A(I,I) + ALPHA*TEMP2
00248                       ELSE
00249                           C(I,J) = BETA*C(I,J) + TEMP1*A(I,I) +
00250      +                             ALPHA*TEMP2
00251                       END IF
00252    90             CONTINUE
00253   100         CONTINUE
00254           END IF
00255       ELSE
00256 *
00257 *        Form  C := alpha*B*A + beta*C.
00258 *
00259           DO 170 J = 1,N
00260               TEMP1 = ALPHA*A(J,J)
00261               IF (BETA.EQ.ZERO) THEN
00262                   DO 110 I = 1,M
00263                       C(I,J) = TEMP1*B(I,J)
00264   110             CONTINUE
00265               ELSE
00266                   DO 120 I = 1,M
00267                       C(I,J) = BETA*C(I,J) + TEMP1*B(I,J)
00268   120             CONTINUE
00269               END IF
00270               DO 140 K = 1,J - 1
00271                   IF (UPPER) THEN
00272                       TEMP1 = ALPHA*A(K,J)
00273                   ELSE
00274                       TEMP1 = ALPHA*A(J,K)
00275                   END IF
00276                   DO 130 I = 1,M
00277                       C(I,J) = C(I,J) + TEMP1*B(I,K)
00278   130             CONTINUE
00279   140         CONTINUE
00280               DO 160 K = J + 1,N
00281                   IF (UPPER) THEN
00282                       TEMP1 = ALPHA*A(J,K)
00283                   ELSE
00284                       TEMP1 = ALPHA*A(K,J)
00285                   END IF
00286                   DO 150 I = 1,M
00287                       C(I,J) = C(I,J) + TEMP1*B(I,K)
00288   150             CONTINUE
00289   160         CONTINUE
00290   170     CONTINUE
00291       END IF
00292 *
00293       RETURN
00294 *
00295 *     End of DSYMM .
00296 *
00297       END
 All Files Functions