git.maemo.org Git - opencv/blob - 3rdparty/lapack/dtrmm.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dtrmm_(char *side, char *uplo, char *transa, char *diag,
   4         integer *m, integer *n, doublereal *alpha, doublereal *a, integer *
   5         lda, doublereal *b, integer *ldb)
   6 {
   7     /* System generated locals */
   8     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
   9
  10     /* Local variables */
  11     integer i__, j, k, info;
  12     doublereal temp;
  13     logical lside;
  14     extern logical lsame_(char *, char *);
  15     integer nrowa;
  16     logical upper;
  17     extern /* Subroutine */ int xerbla_(char *, integer *);
  18     logical nounit;
  19
  20 /*     .. Scalar Arguments .. */
  21 /*     .. */
  22 /*     .. Array Arguments .. */
  23 /*     .. */
  24
  25 /*  Purpose */
  26 /*  ======= */
  27
  28 /*  DTRMM  performs one of the matrix-matrix operations */
  29
  30 /*     B := alpha*op( A )*B,   or   B := alpha*B*op( A ), */
  31
  32 /*  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or */
  33 /*  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of */
  34
  35 /*     op( A ) = A   or   op( A ) = A'. */
  36
  37 /*  Arguments */
  38 /*  ========== */
  39
  40 /*  SIDE   - CHARACTER*1. */
  41 /*           On entry,  SIDE specifies whether  op( A ) multiplies B from */
  42 /*           the left or right as follows: */
  43
  44 /*              SIDE = 'L' or 'l'   B := alpha*op( A )*B. */
  45
  46 /*              SIDE = 'R' or 'r'   B := alpha*B*op( A ). */
  47
  48 /*           Unchanged on exit. */
  49
  50 /*  UPLO   - CHARACTER*1. */
  51 /*           On entry, UPLO specifies whether the matrix A is an upper or */
  52 /*           lower triangular matrix as follows: */
  53
  54 /*              UPLO = 'U' or 'u'   A is an upper triangular matrix. */
  55
  56 /*              UPLO = 'L' or 'l'   A is a lower triangular matrix. */
  57
  58 /*           Unchanged on exit. */
  59
  60 /*  TRANSA - CHARACTER*1. */
  61 /*           On entry, TRANSA specifies the form of op( A ) to be used in */
  62 /*           the matrix multiplication as follows: */
  63
  64 /*              TRANSA = 'N' or 'n'   op( A ) = A. */
  65
  66 /*              TRANSA = 'T' or 't'   op( A ) = A'. */
  67
  68 /*              TRANSA = 'C' or 'c'   op( A ) = A'. */
  69
  70 /*           Unchanged on exit. */
  71
  72 /*  DIAG   - CHARACTER*1. */
  73 /*           On entry, DIAG specifies whether or not A is unit triangular */
  74 /*           as follows: */
  75
  76 /*              DIAG = 'U' or 'u'   A is assumed to be unit triangular. */
  77
  78 /*              DIAG = 'N' or 'n'   A is not assumed to be unit */
  79 /*                                  triangular. */
  80
  81 /*           Unchanged on exit. */
  82
  83 /*  M      - INTEGER. */
  84 /*           On entry, M specifies the number of rows of B. M must be at */
  85 /*           least zero. */
  86 /*           Unchanged on exit. */
  87
  88 /*  N      - INTEGER. */
  89 /*           On entry, N specifies the number of columns of B.  N must be */
  90 /*           at least zero. */
  91 /*           Unchanged on exit. */
  92
  93 /*  ALPHA  - DOUBLE PRECISION. */
  94 /*           On entry,  ALPHA specifies the scalar  alpha. When  alpha is */
  95 /*           zero then  A is not referenced and  B need not be set before */
  96 /*           entry. */
  97 /*           Unchanged on exit. */
  98
  99 /*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m */
 100 /*           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'. */
 101 /*           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k */
 102 /*           upper triangular part of the array  A must contain the upper */
 103 /*           triangular matrix  and the strictly lower triangular part of */
 104 /*           A is not referenced. */
 105 /*           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k */
 106 /*           lower triangular part of the array  A must contain the lower */
 107 /*           triangular matrix  and the strictly upper triangular part of */
 108 /*           A is not referenced. */
 109 /*           Note that when  DIAG = 'U' or 'u',  the diagonal elements of */
 110 /*           A  are not referenced either,  but are assumed to be  unity. */
 111 /*           Unchanged on exit. */
 112
 113 /*  LDA    - INTEGER. */
 114 /*           On entry, LDA specifies the first dimension of A as declared */
 115 /*           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then */
 116 /*           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r' */
 117 /*           then LDA must be at least max( 1, n ). */
 118 /*           Unchanged on exit. */
 119
 120 /*  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ). */
 121 /*           Before entry,  the leading  m by n part of the array  B must */
 122 /*           contain the matrix  B,  and  on exit  is overwritten  by the */
 123 /*           transformed matrix. */
 124
 125 /*  LDB    - INTEGER. */
 126 /*           On entry, LDB specifies the first dimension of B as declared */
 127 /*           in  the  calling  (sub)  program.   LDB  must  be  at  least */
 128 /*           max( 1, m ). */
 129 /*           Unchanged on exit. */
 130
 131
 132 /*  Level 3 Blas routine. */
 133
 134 /*  -- Written on 8-February-1989. */
 135 /*     Jack Dongarra, Argonne National Laboratory. */
 136 /*     Iain Duff, AERE Harwell. */
 137 /*     Jeremy Du Croz, Numerical Algorithms Group Ltd. */
 138 /*     Sven Hammarling, Numerical Algorithms Group Ltd. */
 139
 140
 141 /*     .. External Functions .. */
 142 /*     .. */
 143 /*     .. External Subroutines .. */
 144 /*     .. */
 145 /*     .. Intrinsic Functions .. */
 146 /*     .. */
 147 /*     .. Local Scalars .. */
 148 /*     .. */
 149 /*     .. Parameters .. */
 150 /*     .. */
 151
 152 /*     Test the input parameters. */
 153
 154     /* Parameter adjustments */
 155     a_dim1 = *lda;
 156     a_offset = 1 + a_dim1;
 157     a -= a_offset;
 158     b_dim1 = *ldb;
 159     b_offset = 1 + b_dim1;
 160     b -= b_offset;
 161
 162     /* Function Body */
 163     lside = lsame_(side, "L");
 164     if (lside) {
 165         nrowa = *m;
 166     } else {
 167         nrowa = *n;
 168     }
 169     nounit = lsame_(diag, "N");
 170     upper = lsame_(uplo, "U");
 171
 172     info = 0;
 173     if (! lside && ! lsame_(side, "R")) {
 174         info = 1;
 175     } else if (! upper && ! lsame_(uplo, "L")) {
 176         info = 2;
 177     } else if (! lsame_(transa, "N") && ! lsame_(transa,
 178              "T") && ! lsame_(transa, "C")) {
 179         info = 3;
 180     } else if (! lsame_(diag, "U") && ! lsame_(diag,
 181             "N")) {
 182         info = 4;
 183     } else if (*m < 0) {
 184         info = 5;
 185     } else if (*n < 0) {
 186         info = 6;
 187     } else if (*lda < max(1,nrowa)) {
 188         info = 9;
 189     } else if (*ldb < max(1,*m)) {
 190         info = 11;
 191     }
 192     if (info != 0) {
 193         xerbla_("DTRMM ", &info);
 194         return 0;
 195     }
 196
 197 /*     Quick return if possible. */
 198
 199     if (*n == 0) {
 200         return 0;
 201     }
 202
 203 /*     And when  alpha.eq.zero. */
 204
 205     if (*alpha == 0.) {
 206         i__1 = *n;
 207         for (j = 1; j <= i__1; ++j) {
 208             i__2 = *m;
 209             for (i__ = 1; i__ <= i__2; ++i__) {
 210                 b[i__ + j * b_dim1] = 0.;
 211 /* L10: */
 212             }
 213 /* L20: */
 214         }
 215         return 0;
 216     }
 217
 218 /*     Start the operations. */
 219
 220     if (lside) {
 221         if (lsame_(transa, "N")) {
 222
 223 /*           Form  B := alpha*A*B. */
 224
 225             if (upper) {
 226                 i__1 = *n;
 227                 for (j = 1; j <= i__1; ++j) {
 228                     i__2 = *m;
 229                     for (k = 1; k <= i__2; ++k) {
 230                         if (b[k + j * b_dim1] != 0.) {
 231                             temp = *alpha * b[k + j * b_dim1];
 232                             i__3 = k - 1;
 233                             for (i__ = 1; i__ <= i__3; ++i__) {
 234                                 b[i__ + j * b_dim1] += temp * a[i__ + k *
 235                                         a_dim1];
 236 /* L30: */
 237                             }
 238                             if (nounit) {
 239                                 temp *= a[k + k * a_dim1];
 240                             }
 241                             b[k + j * b_dim1] = temp;
 242                         }
 243 /* L40: */
 244                     }
 245 /* L50: */
 246                 }
 247             } else {
 248                 i__1 = *n;
 249                 for (j = 1; j <= i__1; ++j) {
 250                     for (k = *m; k >= 1; --k) {
 251                         if (b[k + j * b_dim1] != 0.) {
 252                             temp = *alpha * b[k + j * b_dim1];
 253                             b[k + j * b_dim1] = temp;
 254                             if (nounit) {
 255                                 b[k + j * b_dim1] *= a[k + k * a_dim1];
 256                             }
 257                             i__2 = *m;
 258                             for (i__ = k + 1; i__ <= i__2; ++i__) {
 259                                 b[i__ + j * b_dim1] += temp * a[i__ + k *
 260                                         a_dim1];
 261 /* L60: */
 262                             }
 263                         }
 264 /* L70: */
 265                     }
 266 /* L80: */
 267                 }
 268             }
 269         } else {
 270
 271 /*           Form  B := alpha*A'*B. */
 272
 273             if (upper) {
 274                 i__1 = *n;
 275                 for (j = 1; j <= i__1; ++j) {
 276                     for (i__ = *m; i__ >= 1; --i__) {
 277                         temp = b[i__ + j * b_dim1];
 278                         if (nounit) {
 279                             temp *= a[i__ + i__ * a_dim1];
 280                         }
 281                         i__2 = i__ - 1;
 282                         for (k = 1; k <= i__2; ++k) {
 283                             temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
 284 /* L90: */
 285                         }
 286                         b[i__ + j * b_dim1] = *alpha * temp;
 287 /* L100: */
 288                     }
 289 /* L110: */
 290                 }
 291             } else {
 292                 i__1 = *n;
 293                 for (j = 1; j <= i__1; ++j) {
 294                     i__2 = *m;
 295                     for (i__ = 1; i__ <= i__2; ++i__) {
 296                         temp = b[i__ + j * b_dim1];
 297                         if (nounit) {
 298                             temp *= a[i__ + i__ * a_dim1];
 299                         }
 300                         i__3 = *m;
 301                         for (k = i__ + 1; k <= i__3; ++k) {
 302                             temp += a[k + i__ * a_dim1] * b[k + j * b_dim1];
 303 /* L120: */
 304                         }
 305                         b[i__ + j * b_dim1] = *alpha * temp;
 306 /* L130: */
 307                     }
 308 /* L140: */
 309                 }
 310             }
 311         }
 312     } else {
 313         if (lsame_(transa, "N")) {
 314
 315 /*           Form  B := alpha*B*A. */
 316
 317             if (upper) {
 318                 for (j = *n; j >= 1; --j) {
 319                     temp = *alpha;
 320                     if (nounit) {
 321                         temp *= a[j + j * a_dim1];
 322                     }
 323                     i__1 = *m;
 324                     for (i__ = 1; i__ <= i__1; ++i__) {
 325                         b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
 326 /* L150: */
 327                     }
 328                     i__1 = j - 1;
 329                     for (k = 1; k <= i__1; ++k) {
 330                         if (a[k + j * a_dim1] != 0.) {
 331                             temp = *alpha * a[k + j * a_dim1];
 332                             i__2 = *m;
 333                             for (i__ = 1; i__ <= i__2; ++i__) {
 334                                 b[i__ + j * b_dim1] += temp * b[i__ + k *
 335                                         b_dim1];
 336 /* L160: */
 337                             }
 338                         }
 339 /* L170: */
 340                     }
 341 /* L180: */
 342                 }
 343             } else {
 344                 i__1 = *n;
 345                 for (j = 1; j <= i__1; ++j) {
 346                     temp = *alpha;
 347                     if (nounit) {
 348                         temp *= a[j + j * a_dim1];
 349                     }
 350                     i__2 = *m;
 351                     for (i__ = 1; i__ <= i__2; ++i__) {
 352                         b[i__ + j * b_dim1] = temp * b[i__ + j * b_dim1];
 353 /* L190: */
 354                     }
 355                     i__2 = *n;
 356                     for (k = j + 1; k <= i__2; ++k) {
 357                         if (a[k + j * a_dim1] != 0.) {
 358                             temp = *alpha * a[k + j * a_dim1];
 359                             i__3 = *m;
 360                             for (i__ = 1; i__ <= i__3; ++i__) {
 361                                 b[i__ + j * b_dim1] += temp * b[i__ + k *
 362                                         b_dim1];
 363 /* L200: */
 364                             }
 365                         }
 366 /* L210: */
 367                     }
 368 /* L220: */
 369                 }
 370             }
 371         } else {
 372
 373 /*           Form  B := alpha*B*A'. */
 374
 375             if (upper) {
 376                 i__1 = *n;
 377                 for (k = 1; k <= i__1; ++k) {
 378                     i__2 = k - 1;
 379                     for (j = 1; j <= i__2; ++j) {
 380                         if (a[j + k * a_dim1] != 0.) {
 381                             temp = *alpha * a[j + k * a_dim1];
 382                             i__3 = *m;
 383                             for (i__ = 1; i__ <= i__3; ++i__) {
 384                                 b[i__ + j * b_dim1] += temp * b[i__ + k *
 385                                         b_dim1];
 386 /* L230: */
 387                             }
 388                         }
 389 /* L240: */
 390                     }
 391                     temp = *alpha;
 392                     if (nounit) {
 393                         temp *= a[k + k * a_dim1];
 394                     }
 395                     if (temp != 1.) {
 396                         i__2 = *m;
 397                         for (i__ = 1; i__ <= i__2; ++i__) {
 398                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
 399 /* L250: */
 400                         }
 401                     }
 402 /* L260: */
 403                 }
 404             } else {
 405                 for (k = *n; k >= 1; --k) {
 406                     i__1 = *n;
 407                     for (j = k + 1; j <= i__1; ++j) {
 408                         if (a[j + k * a_dim1] != 0.) {
 409                             temp = *alpha * a[j + k * a_dim1];
 410                             i__2 = *m;
 411                             for (i__ = 1; i__ <= i__2; ++i__) {
 412                                 b[i__ + j * b_dim1] += temp * b[i__ + k *
 413                                         b_dim1];
 414 /* L270: */
 415                             }
 416                         }
 417 /* L280: */
 418                     }
 419                     temp = *alpha;
 420                     if (nounit) {
 421                         temp *= a[k + k * a_dim1];
 422                     }
 423                     if (temp != 1.) {
 424                         i__1 = *m;
 425                         for (i__ = 1; i__ <= i__1; ++i__) {
 426                             b[i__ + k * b_dim1] = temp * b[i__ + k * b_dim1];
 427 /* L290: */
 428                         }
 429                     }
 430 /* L300: */
 431                 }
 432             }
 433         }
 434     }
 435
 436     return 0;
 437
 438 /*     End of DTRMM . */
 439
 440 } /* dtrmm_ */