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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7 static real c_b33 = 0.f;
   8 static integer c__0 = 0;
   9
  10 /* Subroutine */ int sgels_(char *trans, integer *m, integer *n, integer *
  11         nrhs, real *a, integer *lda, real *b, integer *ldb, real *work,
  12         integer *lwork, integer *info)
  13 {
  14     /* System generated locals */
  15     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
  16
  17     /* Local variables */
  18     integer i__, j, nb, mn;
  19     real anrm, bnrm;
  20     integer brow;
  21     logical tpsd;
  22     integer iascl, ibscl;
  23     extern logical lsame_(char *, char *);
  24     integer wsize;
  25     real rwork[1];
  26     extern /* Subroutine */ int slabad_(real *, real *);
  27     extern doublereal slamch_(char *), slange_(char *, integer *,
  28             integer *, real *, integer *, real *);
  29     extern /* Subroutine */ int xerbla_(char *, integer *);
  30     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  31             integer *, integer *);
  32     integer scllen;
  33     real bignum;
  34     extern /* Subroutine */ int sgelqf_(integer *, integer *, real *, integer
  35             *, real *, real *, integer *, integer *), slascl_(char *, integer
  36             *, integer *, real *, real *, integer *, integer *, real *,
  37             integer *, integer *), sgeqrf_(integer *, integer *, real
  38             *, integer *, real *, real *, integer *, integer *), slaset_(char
  39             *, integer *, integer *, real *, real *, real *, integer *);
  40     real smlnum;
  41     extern /* Subroutine */ int sormlq_(char *, char *, integer *, integer *,
  42             integer *, real *, integer *, real *, real *, integer *, real *,
  43             integer *, integer *);
  44     logical lquery;
  45     extern /* Subroutine */ int sormqr_(char *, char *, integer *, integer *,
  46             integer *, real *, integer *, real *, real *, integer *, real *,
  47             integer *, integer *), strtrs_(char *, char *,
  48             char *, integer *, integer *, real *, integer *, real *, integer *
  49 , integer *);
  50
  51
  52 /*  -- LAPACK driver routine (version 3.1) -- */
  53 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  54 /*     November 2006 */
  55
  56 /*     .. Scalar Arguments .. */
  57 /*     .. */
  58 /*     .. Array Arguments .. */
  59 /*     .. */
  60
  61 /*  Purpose */
  62 /*  ======= */
  63
  64 /*  SGELS solves overdetermined or underdetermined real linear systems */
  65 /*  involving an M-by-N matrix A, or its transpose, using a QR or LQ */
  66 /*  factorization of A.  It is assumed that A has full rank. */
  67
  68 /*  The following options are provided: */
  69
  70 /*  1. If TRANS = 'N' and m >= n:  find the least squares solution of */
  71 /*     an overdetermined system, i.e., solve the least squares problem */
  72 /*                  minimize || B - A*X ||. */
  73
  74 /*  2. If TRANS = 'N' and m < n:  find the minimum norm solution of */
  75 /*     an underdetermined system A * X = B. */
  76
  77 /*  3. If TRANS = 'T' and m >= n:  find the minimum norm solution of */
  78 /*     an undetermined system A**T * X = B. */
  79
  80 /*  4. If TRANS = 'T' and m < n:  find the least squares solution of */
  81 /*     an overdetermined system, i.e., solve the least squares problem */
  82 /*                  minimize || B - A**T * X ||. */
  83
  84 /*  Several right hand side vectors b and solution vectors x can be */
  85 /*  handled in a single call; they are stored as the columns of the */
  86 /*  M-by-NRHS right hand side matrix B and the N-by-NRHS solution */
  87 /*  matrix X. */
  88
  89 /*  Arguments */
  90 /*  ========= */
  91
  92 /*  TRANS   (input) CHARACTER*1 */
  93 /*          = 'N': the linear system involves A; */
  94 /*          = 'T': the linear system involves A**T. */
  95
  96 /*  M       (input) INTEGER */
  97 /*          The number of rows of the matrix A.  M >= 0. */
  98
  99 /*  N       (input) INTEGER */
 100 /*          The number of columns of the matrix A.  N >= 0. */
 101
 102 /*  NRHS    (input) INTEGER */
 103 /*          The number of right hand sides, i.e., the number of */
 104 /*          columns of the matrices B and X. NRHS >=0. */
 105
 106 /*  A       (input/output) REAL array, dimension (LDA,N) */
 107 /*          On entry, the M-by-N matrix A. */
 108 /*          On exit, */
 109 /*            if M >= N, A is overwritten by details of its QR */
 110 /*                       factorization as returned by SGEQRF; */
 111 /*            if M <  N, A is overwritten by details of its LQ */
 112 /*                       factorization as returned by SGELQF. */
 113
 114 /*  LDA     (input) INTEGER */
 115 /*          The leading dimension of the array A.  LDA >= max(1,M). */
 116
 117 /*  B       (input/output) REAL array, dimension (LDB,NRHS) */
 118 /*          On entry, the matrix B of right hand side vectors, stored */
 119 /*          columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS */
 120 /*          if TRANS = 'T'. */
 121 /*          On exit, if INFO = 0, B is overwritten by the solution */
 122 /*          vectors, stored columnwise: */
 123 /*          if TRANS = 'N' and m >= n, rows 1 to n of B contain the least */
 124 /*          squares solution vectors; the residual sum of squares for the */
 125 /*          solution in each column is given by the sum of squares of */
 126 /*          elements N+1 to M in that column; */
 127 /*          if TRANS = 'N' and m < n, rows 1 to N of B contain the */
 128 /*          minimum norm solution vectors; */
 129 /*          if TRANS = 'T' and m >= n, rows 1 to M of B contain the */
 130 /*          minimum norm solution vectors; */
 131 /*          if TRANS = 'T' and m < n, rows 1 to M of B contain the */
 132 /*          least squares solution vectors; the residual sum of squares */
 133 /*          for the solution in each column is given by the sum of */
 134 /*          squares of elements M+1 to N in that column. */
 135
 136 /*  LDB     (input) INTEGER */
 137 /*          The leading dimension of the array B. LDB >= MAX(1,M,N). */
 138
 139 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
 140 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
 141
 142 /*  LWORK   (input) INTEGER */
 143 /*          The dimension of the array WORK. */
 144 /*          LWORK >= max( 1, MN + max( MN, NRHS ) ). */
 145 /*          For optimal performance, */
 146 /*          LWORK >= max( 1, MN + max( MN, NRHS )*NB ). */
 147 /*          where MN = min(M,N) and NB is the optimum block size. */
 148
 149 /*          If LWORK = -1, then a workspace query is assumed; the routine */
 150 /*          only calculates the optimal size of the WORK array, returns */
 151 /*          this value as the first entry of the WORK array, and no error */
 152 /*          message related to LWORK is issued by XERBLA. */
 153
 154 /*  INFO    (output) INTEGER */
 155 /*          = 0:  successful exit */
 156 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
 157 /*          > 0:  if INFO =  i, the i-th diagonal element of the */
 158 /*                triangular factor of A is zero, so that A does not have */
 159 /*                full rank; the least squares solution could not be */
 160 /*                computed. */
 161
 162 /*  ===================================================================== */
 163
 164 /*     .. Parameters .. */
 165 /*     .. */
 166 /*     .. Local Scalars .. */
 167 /*     .. */
 168 /*     .. Local Arrays .. */
 169 /*     .. */
 170 /*     .. External Functions .. */
 171 /*     .. */
 172 /*     .. External Subroutines .. */
 173 /*     .. */
 174 /*     .. Intrinsic Functions .. */
 175 /*     .. */
 176 /*     .. Executable Statements .. */
 177
 178 /*     Test the input arguments. */
 179
 180     /* Parameter adjustments */
 181     a_dim1 = *lda;
 182     a_offset = 1 + a_dim1;
 183     a -= a_offset;
 184     b_dim1 = *ldb;
 185     b_offset = 1 + b_dim1;
 186     b -= b_offset;
 187     --work;
 188
 189     /* Function Body */
 190     *info = 0;
 191     mn = min(*m,*n);
 192     lquery = *lwork == -1;
 193     if (! (lsame_(trans, "N") || lsame_(trans, "T"))) {
 194         *info = -1;
 195     } else if (*m < 0) {
 196         *info = -2;
 197     } else if (*n < 0) {
 198         *info = -3;
 199     } else if (*nrhs < 0) {
 200         *info = -4;
 201     } else if (*lda < max(1,*m)) {
 202         *info = -6;
 203     } else /* if(complicated condition) */ {
 204 /* Computing MAX */
 205         i__1 = max(1,*m);
 206         if (*ldb < max(i__1,*n)) {
 207             *info = -8;
 208         } else /* if(complicated condition) */ {
 209 /* Computing MAX */
 210             i__1 = 1, i__2 = mn + max(mn,*nrhs);
 211             if (*lwork < max(i__1,i__2) && ! lquery) {
 212                 *info = -10;
 213             }
 214         }
 215     }
 216
 217 /*     Figure out optimal block size */
 218
 219     if (*info == 0 || *info == -10) {
 220
 221         tpsd = TRUE_;
 222         if (lsame_(trans, "N")) {
 223             tpsd = FALSE_;
 224         }
 225
 226         if (*m >= *n) {
 227             nb = ilaenv_(&c__1, "SGEQRF", " ", m, n, &c_n1, &c_n1);
 228             if (tpsd) {
 229 /* Computing MAX */
 230                 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMQR", "LN", m, nrhs, n, &
 231                         c_n1);
 232                 nb = max(i__1,i__2);
 233             } else {
 234 /* Computing MAX */
 235                 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMQR", "LT", m, nrhs, n, &
 236                         c_n1);
 237                 nb = max(i__1,i__2);
 238             }
 239         } else {
 240             nb = ilaenv_(&c__1, "SGELQF", " ", m, n, &c_n1, &c_n1);
 241             if (tpsd) {
 242 /* Computing MAX */
 243                 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMLQ", "LT", n, nrhs, m, &
 244                         c_n1);
 245                 nb = max(i__1,i__2);
 246             } else {
 247 /* Computing MAX */
 248                 i__1 = nb, i__2 = ilaenv_(&c__1, "SORMLQ", "LN", n, nrhs, m, &
 249                         c_n1);
 250                 nb = max(i__1,i__2);
 251             }
 252         }
 253
 254 /* Computing MAX */
 255         i__1 = 1, i__2 = mn + max(mn,*nrhs) * nb;
 256         wsize = max(i__1,i__2);
 257         work[1] = (real) wsize;
 258
 259     }
 260
 261     if (*info != 0) {
 262         i__1 = -(*info);
 263         xerbla_("SGELS ", &i__1);
 264         return 0;
 265     } else if (lquery) {
 266         return 0;
 267     }
 268
 269 /*     Quick return if possible */
 270
 271 /* Computing MIN */
 272     i__1 = min(*m,*n);
 273     if (min(i__1,*nrhs) == 0) {
 274         i__1 = max(*m,*n);
 275         slaset_("Full", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
 276         return 0;
 277     }
 278
 279 /*     Get machine parameters */
 280
 281     smlnum = slamch_("S") / slamch_("P");
 282     bignum = 1.f / smlnum;
 283     slabad_(&smlnum, &bignum);
 284
 285 /*     Scale A, B if max element outside range [SMLNUM,BIGNUM] */
 286
 287     anrm = slange_("M", m, n, &a[a_offset], lda, rwork);
 288     iascl = 0;
 289     if (anrm > 0.f && anrm < smlnum) {
 290
 291 /*        Scale matrix norm up to SMLNUM */
 292
 293         slascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda,
 294                 info);
 295         iascl = 1;
 296     } else if (anrm > bignum) {
 297
 298 /*        Scale matrix norm down to BIGNUM */
 299
 300         slascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda,
 301                 info);
 302         iascl = 2;
 303     } else if (anrm == 0.f) {
 304
 305 /*        Matrix all zero. Return zero solution. */
 306
 307         i__1 = max(*m,*n);
 308         slaset_("F", &i__1, nrhs, &c_b33, &c_b33, &b[b_offset], ldb);
 309         goto L50;
 310     }
 311
 312     brow = *m;
 313     if (tpsd) {
 314         brow = *n;
 315     }
 316     bnrm = slange_("M", &brow, nrhs, &b[b_offset], ldb, rwork);
 317     ibscl = 0;
 318     if (bnrm > 0.f && bnrm < smlnum) {
 319
 320 /*        Scale matrix norm up to SMLNUM */
 321
 322         slascl_("G", &c__0, &c__0, &bnrm, &smlnum, &brow, nrhs, &b[b_offset],
 323                 ldb, info);
 324         ibscl = 1;
 325     } else if (bnrm > bignum) {
 326
 327 /*        Scale matrix norm down to BIGNUM */
 328
 329         slascl_("G", &c__0, &c__0, &bnrm, &bignum, &brow, nrhs, &b[b_offset],
 330                 ldb, info);
 331         ibscl = 2;
 332     }
 333
 334     if (*m >= *n) {
 335
 336 /*        compute QR factorization of A */
 337
 338         i__1 = *lwork - mn;
 339         sgeqrf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
 340                 ;
 341
 342 /*        workspace at least N, optimally N*NB */
 343
 344         if (! tpsd) {
 345
 346 /*           Least-Squares Problem min || A * X - B || */
 347
 348 /*           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */
 349
 350             i__1 = *lwork - mn;
 351             sormqr_("Left", "Transpose", m, nrhs, n, &a[a_offset], lda, &work[
 352                     1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
 353
 354 /*           workspace at least NRHS, optimally NRHS*NB */
 355
 356 /*           B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS) */
 357
 358             strtrs_("Upper", "No transpose", "Non-unit", n, nrhs, &a[a_offset]
 359 , lda, &b[b_offset], ldb, info);
 360
 361             if (*info > 0) {
 362                 return 0;
 363             }
 364
 365             scllen = *n;
 366
 367         } else {
 368
 369 /*           Overdetermined system of equations A' * X = B */
 370
 371 /*           B(1:N,1:NRHS) := inv(R') * B(1:N,1:NRHS) */
 372
 373             strtrs_("Upper", "Transpose", "Non-unit", n, nrhs, &a[a_offset],
 374                     lda, &b[b_offset], ldb, info);
 375
 376             if (*info > 0) {
 377                 return 0;
 378             }
 379
 380 /*           B(N+1:M,1:NRHS) = ZERO */
 381
 382             i__1 = *nrhs;
 383             for (j = 1; j <= i__1; ++j) {
 384                 i__2 = *m;
 385                 for (i__ = *n + 1; i__ <= i__2; ++i__) {
 386                     b[i__ + j * b_dim1] = 0.f;
 387 /* L10: */
 388                 }
 389 /* L20: */
 390             }
 391
 392 /*           B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS) */
 393
 394             i__1 = *lwork - mn;
 395             sormqr_("Left", "No transpose", m, nrhs, n, &a[a_offset], lda, &
 396                     work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
 397
 398 /*           workspace at least NRHS, optimally NRHS*NB */
 399
 400             scllen = *m;
 401
 402         }
 403
 404     } else {
 405
 406 /*        Compute LQ factorization of A */
 407
 408         i__1 = *lwork - mn;
 409         sgelqf_(m, n, &a[a_offset], lda, &work[1], &work[mn + 1], &i__1, info)
 410                 ;
 411
 412 /*        workspace at least M, optimally M*NB. */
 413
 414         if (! tpsd) {
 415
 416 /*           underdetermined system of equations A * X = B */
 417
 418 /*           B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS) */
 419
 420             strtrs_("Lower", "No transpose", "Non-unit", m, nrhs, &a[a_offset]
 421 , lda, &b[b_offset], ldb, info);
 422
 423             if (*info > 0) {
 424                 return 0;
 425             }
 426
 427 /*           B(M+1:N,1:NRHS) = 0 */
 428
 429             i__1 = *nrhs;
 430             for (j = 1; j <= i__1; ++j) {
 431                 i__2 = *n;
 432                 for (i__ = *m + 1; i__ <= i__2; ++i__) {
 433                     b[i__ + j * b_dim1] = 0.f;
 434 /* L30: */
 435                 }
 436 /* L40: */
 437             }
 438
 439 /*           B(1:N,1:NRHS) := Q(1:N,:)' * B(1:M,1:NRHS) */
 440
 441             i__1 = *lwork - mn;
 442             sormlq_("Left", "Transpose", n, nrhs, m, &a[a_offset], lda, &work[
 443                     1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
 444
 445 /*           workspace at least NRHS, optimally NRHS*NB */
 446
 447             scllen = *n;
 448
 449         } else {
 450
 451 /*           overdetermined system min || A' * X - B || */
 452
 453 /*           B(1:N,1:NRHS) := Q * B(1:N,1:NRHS) */
 454
 455             i__1 = *lwork - mn;
 456             sormlq_("Left", "No transpose", n, nrhs, m, &a[a_offset], lda, &
 457                     work[1], &b[b_offset], ldb, &work[mn + 1], &i__1, info);
 458
 459 /*           workspace at least NRHS, optimally NRHS*NB */
 460
 461 /*           B(1:M,1:NRHS) := inv(L') * B(1:M,1:NRHS) */
 462
 463             strtrs_("Lower", "Transpose", "Non-unit", m, nrhs, &a[a_offset],
 464                     lda, &b[b_offset], ldb, info);
 465
 466             if (*info > 0) {
 467                 return 0;
 468             }
 469
 470             scllen = *m;
 471
 472         }
 473
 474     }
 475
 476 /*     Undo scaling */
 477
 478     if (iascl == 1) {
 479         slascl_("G", &c__0, &c__0, &anrm, &smlnum, &scllen, nrhs, &b[b_offset]
 480 , ldb, info);
 481     } else if (iascl == 2) {
 482         slascl_("G", &c__0, &c__0, &anrm, &bignum, &scllen, nrhs, &b[b_offset]
 483 , ldb, info);
 484     }
 485     if (ibscl == 1) {
 486         slascl_("G", &c__0, &c__0, &smlnum, &bnrm, &scllen, nrhs, &b[b_offset]
 487 , ldb, info);
 488     } else if (ibscl == 2) {
 489         slascl_("G", &c__0, &c__0, &bignum, &bnrm, &scllen, nrhs, &b[b_offset]
 490 , ldb, info);
 491     }
 492
 493 L50:
 494     work[1] = (real) wsize;
 495
 496     return 0;
 497
 498 /*     End of SGELS */
 499
 500 } /* sgels_ */