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

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