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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dsytrs_(char *uplo, integer *n, integer *nrhs,
   4         doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *
   5         ldb, integer *info)
   6 {
   7 /*  -- LAPACK routine (version 3.0) --
   8        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
   9        Courant Institute, Argonne National Lab, and Rice University
  10        March 31, 1993
  11
  12
  13     Purpose
  14     =======
  15
  16     DSYTRS solves a system of linear equations A*X = B with a real
  17     symmetric matrix A using the factorization A = U*D*U**T or
  18     A = L*D*L**T computed by DSYTRF.
  19
  20     Arguments
  21     =========
  22
  23     UPLO    (input) CHARACTER*1
  24             Specifies whether the details of the factorization are stored
  25             as an upper or lower triangular matrix.
  26             = 'U':  Upper triangular, form is A = U*D*U**T;
  27             = 'L':  Lower triangular, form is A = L*D*L**T.
  28
  29     N       (input) INTEGER
  30             The order of the matrix A.  N >= 0.
  31
  32     NRHS    (input) INTEGER
  33             The number of right hand sides, i.e., the number of columns
  34             of the matrix B.  NRHS >= 0.
  35
  36     A       (input) DOUBLE PRECISION array, dimension (LDA,N)
  37             The block diagonal matrix D and the multipliers used to
  38             obtain the factor U or L as computed by DSYTRF.
  39
  40     LDA     (input) INTEGER
  41             The leading dimension of the array A.  LDA >= max(1,N).
  42
  43     IPIV    (input) INTEGER array, dimension (N)
  44             Details of the interchanges and the block structure of D
  45             as determined by DSYTRF.
  46
  47     B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
  48             On entry, the right hand side matrix B.
  49             On exit, the solution matrix X.
  50
  51     LDB     (input) INTEGER
  52             The leading dimension of the array B.  LDB >= max(1,N).
  53
  54     INFO    (output) INTEGER
  55             = 0:  successful exit
  56             < 0:  if INFO = -i, the i-th argument had an illegal value
  57
  58     =====================================================================
  59
  60
  61        Parameter adjustments */
  62     /* Table of constant values */
  63     static doublereal c_b7 = -1.;
  64     static integer c__1 = 1;
  65     static doublereal c_b19 = 1.;
  66
  67     /* System generated locals */
  68     integer a_dim1, a_offset, b_dim1, b_offset, i__1;
  69     doublereal d__1;
  70     /* Local variables */
  71     extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
  72             doublereal *, integer *, doublereal *, integer *, doublereal *,
  73             integer *);
  74     static doublereal akm1k;
  75     static integer j, k;
  76     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  77             integer *);
  78     extern logical lsame_(char *, char *);
  79     static doublereal denom;
  80     extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
  81             doublereal *, doublereal *, integer *, doublereal *, integer *,
  82             doublereal *, doublereal *, integer *), dswap_(integer *,
  83             doublereal *, integer *, doublereal *, integer *);
  84     static logical upper;
  85     static doublereal ak, bk;
  86     static integer kp;
  87     extern /* Subroutine */ int xerbla_(char *, integer *);
  88     static doublereal akm1, bkm1;
  89 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
  90 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
  91
  92
  93     a_dim1 = *lda;
  94     a_offset = 1 + a_dim1 * 1;
  95     a -= a_offset;
  96     --ipiv;
  97     b_dim1 = *ldb;
  98     b_offset = 1 + b_dim1 * 1;
  99     b -= b_offset;
 100
 101     /* Function Body */
 102     *info = 0;
 103     upper = lsame_(uplo, "U");
 104     if (! upper && ! lsame_(uplo, "L")) {
 105         *info = -1;
 106     } else if (*n < 0) {
 107         *info = -2;
 108     } else if (*nrhs < 0) {
 109         *info = -3;
 110     } else if (*lda < max(1,*n)) {
 111         *info = -5;
 112     } else if (*ldb < max(1,*n)) {
 113         *info = -8;
 114     }
 115     if (*info != 0) {
 116         i__1 = -(*info);
 117         xerbla_("DSYTRS", &i__1);
 118         return 0;
 119     }
 120
 121 /*     Quick return if possible */
 122
 123     if (*n == 0 || *nrhs == 0) {
 124         return 0;
 125     }
 126
 127     if (upper) {
 128
 129 /*        Solve A*X = B, where A = U*D*U'.
 130
 131           First solve U*D*X = B, overwriting B with X.
 132
 133           K is the main loop index, decreasing from N to 1 in steps of
 134           1 or 2, depending on the size of the diagonal blocks. */
 135
 136         k = *n;
 137 L10:
 138
 139 /*        If K < 1, exit from loop. */
 140
 141         if (k < 1) {
 142             goto L30;
 143         }
 144
 145         if (ipiv[k] > 0) {
 146
 147 /*           1 x 1 diagonal block
 148
 149              Interchange rows K and IPIV(K). */
 150
 151             kp = ipiv[k];
 152             if (kp != k) {
 153                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 154             }
 155
 156 /*           Multiply by inv(U(K)), where U(K) is the transformation
 157              stored in column K of A. */
 158
 159             i__1 = k - 1;
 160             dger_(&i__1, nrhs, &c_b7, &a_ref(1, k), &c__1, &b_ref(k, 1), ldb,
 161                     &b_ref(1, 1), ldb);
 162
 163 /*           Multiply by the inverse of the diagonal block. */
 164
 165             d__1 = 1. / a_ref(k, k);
 166             dscal_(nrhs, &d__1, &b_ref(k, 1), ldb);
 167             --k;
 168         } else {
 169
 170 /*           2 x 2 diagonal block
 171
 172              Interchange rows K-1 and -IPIV(K). */
 173
 174             kp = -ipiv[k];
 175             if (kp != k - 1) {
 176                 dswap_(nrhs, &b_ref(k - 1, 1), ldb, &b_ref(kp, 1), ldb);
 177             }
 178
 179 /*           Multiply by inv(U(K)), where U(K) is the transformation
 180              stored in columns K-1 and K of A. */
 181
 182             i__1 = k - 2;
 183             dger_(&i__1, nrhs, &c_b7, &a_ref(1, k), &c__1, &b_ref(k, 1), ldb,
 184                     &b_ref(1, 1), ldb);
 185             i__1 = k - 2;
 186             dger_(&i__1, nrhs, &c_b7, &a_ref(1, k - 1), &c__1, &b_ref(k - 1,
 187                     1), ldb, &b_ref(1, 1), ldb);
 188
 189 /*           Multiply by the inverse of the diagonal block. */
 190
 191             akm1k = a_ref(k - 1, k);
 192             akm1 = a_ref(k - 1, k - 1) / akm1k;
 193             ak = a_ref(k, k) / akm1k;
 194             denom = akm1 * ak - 1.;
 195             i__1 = *nrhs;
 196             for (j = 1; j <= i__1; ++j) {
 197                 bkm1 = b_ref(k - 1, j) / akm1k;
 198                 bk = b_ref(k, j) / akm1k;
 199                 b_ref(k - 1, j) = (ak * bkm1 - bk) / denom;
 200                 b_ref(k, j) = (akm1 * bk - bkm1) / denom;
 201 /* L20: */
 202             }
 203             k += -2;
 204         }
 205
 206         goto L10;
 207 L30:
 208
 209 /*        Next solve U'*X = B, overwriting B with X.
 210
 211           K is the main loop index, increasing from 1 to N in steps of
 212           1 or 2, depending on the size of the diagonal blocks. */
 213
 214         k = 1;
 215 L40:
 216
 217 /*        If K > N, exit from loop. */
 218
 219         if (k > *n) {
 220             goto L50;
 221         }
 222
 223         if (ipiv[k] > 0) {
 224
 225 /*           1 x 1 diagonal block
 226
 227              Multiply by inv(U'(K)), where U(K) is the transformation
 228              stored in column K of A. */
 229
 230             i__1 = k - 1;
 231             dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &a_ref(
 232                     1, k), &c__1, &c_b19, &b_ref(k, 1), ldb);
 233
 234 /*           Interchange rows K and IPIV(K). */
 235
 236             kp = ipiv[k];
 237             if (kp != k) {
 238                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 239             }
 240             ++k;
 241         } else {
 242
 243 /*           2 x 2 diagonal block
 244
 245              Multiply by inv(U'(K+1)), where U(K+1) is the transformation
 246              stored in columns K and K+1 of A. */
 247
 248             i__1 = k - 1;
 249             dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &a_ref(
 250                     1, k), &c__1, &c_b19, &b_ref(k, 1), ldb);
 251             i__1 = k - 1;
 252             dgemv_("Transpose", &i__1, nrhs, &c_b7, &b[b_offset], ldb, &a_ref(
 253                     1, k + 1), &c__1, &c_b19, &b_ref(k + 1, 1), ldb);
 254
 255 /*           Interchange rows K and -IPIV(K). */
 256
 257             kp = -ipiv[k];
 258             if (kp != k) {
 259                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 260             }
 261             k += 2;
 262         }
 263
 264         goto L40;
 265 L50:
 266
 267         ;
 268     } else {
 269
 270 /*        Solve A*X = B, where A = L*D*L'.
 271
 272           First solve L*D*X = B, overwriting B with X.
 273
 274           K is the main loop index, increasing from 1 to N in steps of
 275           1 or 2, depending on the size of the diagonal blocks. */
 276
 277         k = 1;
 278 L60:
 279
 280 /*        If K > N, exit from loop. */
 281
 282         if (k > *n) {
 283             goto L80;
 284         }
 285
 286         if (ipiv[k] > 0) {
 287
 288 /*           1 x 1 diagonal block
 289
 290              Interchange rows K and IPIV(K). */
 291
 292             kp = ipiv[k];
 293             if (kp != k) {
 294                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 295             }
 296
 297 /*           Multiply by inv(L(K)), where L(K) is the transformation
 298              stored in column K of A. */
 299
 300             if (k < *n) {
 301                 i__1 = *n - k;
 302                 dger_(&i__1, nrhs, &c_b7, &a_ref(k + 1, k), &c__1, &b_ref(k,
 303                         1), ldb, &b_ref(k + 1, 1), ldb);
 304             }
 305
 306 /*           Multiply by the inverse of the diagonal block. */
 307
 308             d__1 = 1. / a_ref(k, k);
 309             dscal_(nrhs, &d__1, &b_ref(k, 1), ldb);
 310             ++k;
 311         } else {
 312
 313 /*           2 x 2 diagonal block
 314
 315              Interchange rows K+1 and -IPIV(K). */
 316
 317             kp = -ipiv[k];
 318             if (kp != k + 1) {
 319                 dswap_(nrhs, &b_ref(k + 1, 1), ldb, &b_ref(kp, 1), ldb);
 320             }
 321
 322 /*           Multiply by inv(L(K)), where L(K) is the transformation
 323              stored in columns K and K+1 of A. */
 324
 325             if (k < *n - 1) {
 326                 i__1 = *n - k - 1;
 327                 dger_(&i__1, nrhs, &c_b7, &a_ref(k + 2, k), &c__1, &b_ref(k,
 328                         1), ldb, &b_ref(k + 2, 1), ldb);
 329                 i__1 = *n - k - 1;
 330                 dger_(&i__1, nrhs, &c_b7, &a_ref(k + 2, k + 1), &c__1, &b_ref(
 331                         k + 1, 1), ldb, &b_ref(k + 2, 1), ldb);
 332             }
 333
 334 /*           Multiply by the inverse of the diagonal block. */
 335
 336             akm1k = a_ref(k + 1, k);
 337             akm1 = a_ref(k, k) / akm1k;
 338             ak = a_ref(k + 1, k + 1) / akm1k;
 339             denom = akm1 * ak - 1.;
 340             i__1 = *nrhs;
 341             for (j = 1; j <= i__1; ++j) {
 342                 bkm1 = b_ref(k, j) / akm1k;
 343                 bk = b_ref(k + 1, j) / akm1k;
 344                 b_ref(k, j) = (ak * bkm1 - bk) / denom;
 345                 b_ref(k + 1, j) = (akm1 * bk - bkm1) / denom;
 346 /* L70: */
 347             }
 348             k += 2;
 349         }
 350
 351         goto L60;
 352 L80:
 353
 354 /*        Next solve L'*X = B, overwriting B with X.
 355
 356           K is the main loop index, decreasing from N to 1 in steps of
 357           1 or 2, depending on the size of the diagonal blocks. */
 358
 359         k = *n;
 360 L90:
 361
 362 /*        If K < 1, exit from loop. */
 363
 364         if (k < 1) {
 365             goto L100;
 366         }
 367
 368         if (ipiv[k] > 0) {
 369
 370 /*           1 x 1 diagonal block
 371
 372              Multiply by inv(L'(K)), where L(K) is the transformation
 373              stored in column K of A. */
 374
 375             if (k < *n) {
 376                 i__1 = *n - k;
 377                 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b_ref(k + 1, 1), ldb,
 378                          &a_ref(k + 1, k), &c__1, &c_b19, &b_ref(k, 1), ldb);
 379             }
 380
 381 /*           Interchange rows K and IPIV(K). */
 382
 383             kp = ipiv[k];
 384             if (kp != k) {
 385                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 386             }
 387             --k;
 388         } else {
 389
 390 /*           2 x 2 diagonal block
 391
 392              Multiply by inv(L'(K-1)), where L(K-1) is the transformation
 393              stored in columns K-1 and K of A. */
 394
 395             if (k < *n) {
 396                 i__1 = *n - k;
 397                 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b_ref(k + 1, 1), ldb,
 398                          &a_ref(k + 1, k), &c__1, &c_b19, &b_ref(k, 1), ldb);
 399                 i__1 = *n - k;
 400                 dgemv_("Transpose", &i__1, nrhs, &c_b7, &b_ref(k + 1, 1), ldb,
 401                          &a_ref(k + 1, k - 1), &c__1, &c_b19, &b_ref(k - 1, 1)
 402                         , ldb);
 403             }
 404
 405 /*           Interchange rows K and -IPIV(K). */
 406
 407             kp = -ipiv[k];
 408             if (kp != k) {
 409                 dswap_(nrhs, &b_ref(k, 1), ldb, &b_ref(kp, 1), ldb);
 410             }
 411             k += -2;
 412         }
 413
 414         goto L90;
 415 L100:
 416         ;
 417     }
 418
 419     return 0;
 420
 421 /*     End of DSYTRS */
 422
 423 } /* dsytrs_ */
 424
 425 #undef b_ref
 426 #undef a_ref
 427
 428