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