sgeqrf.c
Go to the documentation of this file.
00001 /* sgeqrf.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static integer c_n1 = -1;
00020 static integer c__3 = 3;
00021 static integer c__2 = 2;
00022 
00023 /* Subroutine */ int sgeqrf_(integer *m, integer *n, real *a, integer *lda, 
00024         real *tau, real *work, integer *lwork, integer *info)
00025 {
00026     /* System generated locals */
00027     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5, i__6;
00028     real r__1;
00029 
00030     /* Local variables */
00031     integer i__, j, k, ib, nb, nt, nx, iws;
00032     extern doublereal sceil_(real *);
00033     integer nbmin, iinfo;
00034     extern /* Subroutine */ int sgeqr2_(integer *, integer *, real *, integer 
00035             *, real *, real *, integer *), slarfb_(char *, char *, char *, 
00036             char *, integer *, integer *, integer *, real *, integer *, real *
00037 , integer *, real *, integer *, real *, integer *), xerbla_(char *, integer *);
00038     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00039             integer *, integer *);
00040     extern /* Subroutine */ int slarft_(char *, char *, integer *, integer *, 
00041             real *, integer *, real *, real *, integer *);
00042     integer lbwork, llwork, lwkopt;
00043     logical lquery;
00044 
00045 
00046 /*  -- LAPACK routine (version 3.1) -- */
00047 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00048 /*     March 2008 */
00049 
00050 /*     .. Scalar Arguments .. */
00051 /*     .. */
00052 /*     .. Array Arguments .. */
00053 /*     .. */
00054 
00055 /*  Purpose */
00056 /*  ======= */
00057 
00058 /*  SGEQRF computes a QR factorization of a real M-by-N matrix A: */
00059 /*  A = Q * R. */
00060 
00061 /*  This is the left-looking Level 3 BLAS version of the algorithm. */
00062 
00063 /*  Arguments */
00064 /*  ========= */
00065 
00066 /*  M       (input) INTEGER */
00067 /*          The number of rows of the matrix A.  M >= 0. */
00068 
00069 /*  N       (input) INTEGER */
00070 /*          The number of columns of the matrix A.  N >= 0. */
00071 
00072 /*  A       (input/output) REAL array, dimension (LDA,N) */
00073 /*          On entry, the M-by-N matrix A. */
00074 /*          On exit, the elements on and above the diagonal of the array */
00075 /*          contain the min(M,N)-by-N upper trapezoidal matrix R (R is */
00076 /*          upper triangular if m >= n); the elements below the diagonal, */
00077 /*          with the array TAU, represent the orthogonal matrix Q as a */
00078 /*          product of min(m,n) elementary reflectors (see Further */
00079 /*          Details). */
00080 
00081 /*  LDA     (input) INTEGER */
00082 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00083 
00084 /*  TAU     (output) REAL array, dimension (min(M,N)) */
00085 /*          The scalar factors of the elementary reflectors (see Further */
00086 /*          Details). */
00087 
00088 /*  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
00089 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00090 
00091 /*  LWORK   (input) INTEGER */
00092 
00093 /*          The dimension of the array WORK. The dimension can be divided into three parts. */
00094 
00095 /*          1) The part for the triangular factor T. If the very last T is not bigger */
00096 /*             than any of the rest, then this part is NB x ceiling(K/NB), otherwise, */
00097 /*             NB x (K-NT), where K = min(M,N) and NT is the dimension of the very last T */
00098 
00099 /*          2) The part for the very last T when T is bigger than any of the rest T. */
00100 /*             The size of this part is NT x NT, where NT = K - ceiling ((K-NX)/NB) x NB, */
00101 /*             where K = min(M,N), NX is calculated by */
00102 /*                   NX = MAX( 0, ILAENV( 3, 'SGEQRF', ' ', M, N, -1, -1 ) ) */
00103 
00104 /*          3) The part for dlarfb is of size max((N-M)*K, (N-M)*NB, K*NB, NB*NB) */
00105 
00106 /*          So LWORK = part1 + part2 + part3 */
00107 
00108 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00109 /*          only calculates the optimal size of the WORK array, returns */
00110 /*          this value as the first entry of the WORK array, and no error */
00111 /*          message related to LWORK is issued by XERBLA. */
00112 
00113 /*  INFO    (output) INTEGER */
00114 /*          = 0:  successful exit */
00115 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00116 
00117 /*  Further Details */
00118 /*  =============== */
00119 
00120 /*  The matrix Q is represented as a product of elementary reflectors */
00121 
00122 /*     Q = H(1) H(2) . . . H(k), where k = min(m,n). */
00123 
00124 /*  Each H(i) has the form */
00125 
00126 /*     H(i) = I - tau * v * v' */
00127 
00128 /*  where tau is a real scalar, and v is a real vector with */
00129 /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i), */
00130 /*  and tau in TAU(i). */
00131 
00132 /*  ===================================================================== */
00133 
00134 /*     .. Local Scalars .. */
00135 /*     .. */
00136 /*     .. External Subroutines .. */
00137 /*     .. */
00138 /*     .. Intrinsic Functions .. */
00139 /*     .. */
00140 /*     .. External Functions .. */
00141 /*     .. */
00142 /*     .. Executable Statements .. */
00143     /* Parameter adjustments */
00144     a_dim1 = *lda;
00145     a_offset = 1 + a_dim1;
00146     a -= a_offset;
00147     --tau;
00148     --work;
00149 
00150     /* Function Body */
00151     *info = 0;
00152     nbmin = 2;
00153     nx = 0;
00154     iws = *n;
00155     k = min(*m,*n);
00156     nb = ilaenv_(&c__1, "SGEQRF", " ", m, n, &c_n1, &c_n1);
00157     if (nb > 1 && nb < k) {
00158 
00159 /*        Determine when to cross over from blocked to unblocked code. */
00160 
00161 /* Computing MAX */
00162         i__1 = 0, i__2 = ilaenv_(&c__3, "SGEQRF", " ", m, n, &c_n1, &c_n1);
00163         nx = max(i__1,i__2);
00164     }
00165 
00166 /*     Get NT, the size of the very last T, which is the left-over from in-between K-NX and K to K, eg.: */
00167 
00168 /*            NB=3     2NB=6       K=10 */
00169 /*            |        |           | */
00170 /*      1--2--3--4--5--6--7--8--9--10 */
00171 /*                  |     \________/ */
00172 /*               K-NX=5      NT=4 */
00173 
00174 /*     So here 4 x 4 is the last T stored in the workspace */
00175 
00176     r__1 = (real) (k - nx) / (real) nb;
00177     nt = k - sceil_(&r__1) * nb;
00178 
00179 /*     optimal workspace = space for dlarfb + space for normal T's + space for the last T */
00180 
00181 /* Computing MAX */
00182 /* Computing MAX */
00183     i__3 = (*n - *m) * k, i__4 = (*n - *m) * nb;
00184 /* Computing MAX */
00185     i__5 = k * nb, i__6 = nb * nb;
00186     i__1 = max(i__3,i__4), i__2 = max(i__5,i__6);
00187     llwork = max(i__1,i__2);
00188     r__1 = (real) llwork / (real) nb;
00189     llwork = sceil_(&r__1);
00190     if (nt > nb) {
00191         lbwork = k - nt;
00192 
00193 /*         Optimal workspace for dlarfb = MAX(1,N)*NT */
00194 
00195         lwkopt = (lbwork + llwork) * nb;
00196         work[1] = (real) (lwkopt + nt * nt);
00197     } else {
00198         r__1 = (real) k / (real) nb;
00199         lbwork = sceil_(&r__1) * nb;
00200         lwkopt = (lbwork + llwork - nb) * nb;
00201         work[1] = (real) lwkopt;
00202     }
00203 
00204 /*     Test the input arguments */
00205 
00206     lquery = *lwork == -1;
00207     if (*m < 0) {
00208         *info = -1;
00209     } else if (*n < 0) {
00210         *info = -2;
00211     } else if (*lda < max(1,*m)) {
00212         *info = -4;
00213     } else if (*lwork < max(1,*n) && ! lquery) {
00214         *info = -7;
00215     }
00216     if (*info != 0) {
00217         i__1 = -(*info);
00218         xerbla_("SGEQRF", &i__1);
00219         return 0;
00220     } else if (lquery) {
00221         return 0;
00222     }
00223 
00224 /*     Quick return if possible */
00225 
00226     if (k == 0) {
00227         work[1] = 1.f;
00228         return 0;
00229     }
00230 
00231     if (nb > 1 && nb < k) {
00232         if (nx < k) {
00233 
00234 /*           Determine if workspace is large enough for blocked code. */
00235 
00236             if (nt <= nb) {
00237                 iws = (lbwork + llwork - nb) * nb;
00238             } else {
00239                 iws = (lbwork + llwork) * nb + nt * nt;
00240             }
00241             if (*lwork < iws) {
00242 
00243 /*              Not enough workspace to use optimal NB:  reduce NB and */
00244 /*              determine the minimum value of NB. */
00245 
00246                 if (nt <= nb) {
00247                     nb = *lwork / (llwork + (lbwork - nb));
00248                 } else {
00249                     nb = (*lwork - nt * nt) / (lbwork + llwork);
00250                 }
00251 /* Computing MAX */
00252                 i__1 = 2, i__2 = ilaenv_(&c__2, "SGEQRF", " ", m, n, &c_n1, &
00253                         c_n1);
00254                 nbmin = max(i__1,i__2);
00255             }
00256         }
00257     }
00258 
00259     if (nb >= nbmin && nb < k && nx < k) {
00260 
00261 /*        Use blocked code initially */
00262 
00263         i__1 = k - nx;
00264         i__2 = nb;
00265         for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
00266 /* Computing MIN */
00267             i__3 = k - i__ + 1;
00268             ib = min(i__3,nb);
00269 
00270 /*           Update the current column using old T's */
00271 
00272             i__3 = i__ - nb;
00273             i__4 = nb;
00274             for (j = 1; i__4 < 0 ? j >= i__3 : j <= i__3; j += i__4) {
00275 
00276 /*              Apply H' to A(J:M,I:I+IB-1) from the left */
00277 
00278                 i__5 = *m - j + 1;
00279                 slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__5, &
00280                         ib, &nb, &a[j + j * a_dim1], lda, &work[j], &lbwork, &
00281                         a[j + i__ * a_dim1], lda, &work[lbwork * nb + nt * nt 
00282                         + 1], &ib);
00283 /* L20: */
00284             }
00285 
00286 /*           Compute the QR factorization of the current block */
00287 /*           A(I:M,I:I+IB-1) */
00288 
00289             i__4 = *m - i__ + 1;
00290             sgeqr2_(&i__4, &ib, &a[i__ + i__ * a_dim1], lda, &tau[i__], &work[
00291                     lbwork * nb + nt * nt + 1], &iinfo);
00292             if (i__ + ib <= *n) {
00293 
00294 /*              Form the triangular factor of the block reflector */
00295 /*              H = H(i) H(i+1) . . . H(i+ib-1) */
00296 
00297                 i__4 = *m - i__ + 1;
00298                 slarft_("Forward", "Columnwise", &i__4, &ib, &a[i__ + i__ * 
00299                         a_dim1], lda, &tau[i__], &work[i__], &lbwork);
00300 
00301             }
00302 /* L10: */
00303         }
00304     } else {
00305         i__ = 1;
00306     }
00307 
00308 /*     Use unblocked code to factor the last or only block. */
00309 
00310     if (i__ <= k) {
00311         if (i__ != 1) {
00312             i__2 = i__ - nb;
00313             i__1 = nb;
00314             for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
00315 
00316 /*                Apply H' to A(J:M,I:K) from the left */
00317 
00318                 i__4 = *m - j + 1;
00319                 i__3 = k - i__ + 1;
00320                 i__5 = k - i__ + 1;
00321                 slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__4, &
00322                         i__3, &nb, &a[j + j * a_dim1], lda, &work[j], &lbwork, 
00323                          &a[j + i__ * a_dim1], lda, &work[lbwork * nb + nt * 
00324                         nt + 1], &i__5);
00325 /* L30: */
00326             }
00327             i__1 = *m - i__ + 1;
00328             i__2 = k - i__ + 1;
00329             sgeqr2_(&i__1, &i__2, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
00330                     work[lbwork * nb + nt * nt + 1], &iinfo);
00331         } else {
00332 
00333 /*        Use unblocked code to factor the last or only block. */
00334 
00335             i__1 = *m - i__ + 1;
00336             i__2 = *n - i__ + 1;
00337             sgeqr2_(&i__1, &i__2, &a[i__ + i__ * a_dim1], lda, &tau[i__], &
00338                     work[1], &iinfo);
00339         }
00340     }
00341 
00342 /*     Apply update to the column M+1:N when N > M */
00343 
00344     if (*m < *n && i__ != 1) {
00345 
00346 /*         Form the last triangular factor of the block reflector */
00347 /*         H = H(i) H(i+1) . . . H(i+ib-1) */
00348 
00349         if (nt <= nb) {
00350             i__1 = *m - i__ + 1;
00351             i__2 = k - i__ + 1;
00352             slarft_("Forward", "Columnwise", &i__1, &i__2, &a[i__ + i__ * 
00353                     a_dim1], lda, &tau[i__], &work[i__], &lbwork);
00354         } else {
00355             i__1 = *m - i__ + 1;
00356             i__2 = k - i__ + 1;
00357             slarft_("Forward", "Columnwise", &i__1, &i__2, &a[i__ + i__ * 
00358                     a_dim1], lda, &tau[i__], &work[lbwork * nb + 1], &nt);
00359         }
00360 
00361 /*         Apply H' to A(1:M,M+1:N) from the left */
00362 
00363         i__1 = k - nx;
00364         i__2 = nb;
00365         for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
00366 /* Computing MIN */
00367             i__4 = k - j + 1;
00368             ib = min(i__4,nb);
00369             i__4 = *m - j + 1;
00370             i__3 = *n - *m;
00371             i__5 = *n - *m;
00372             slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__4, &
00373                     i__3, &ib, &a[j + j * a_dim1], lda, &work[j], &lbwork, &a[
00374                     j + (*m + 1) * a_dim1], lda, &work[lbwork * nb + nt * nt 
00375                     + 1], &i__5);
00376 /* L40: */
00377         }
00378         if (nt <= nb) {
00379             i__2 = *m - j + 1;
00380             i__1 = *n - *m;
00381             i__4 = k - j + 1;
00382             i__3 = *n - *m;
00383             slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__2, &
00384                     i__1, &i__4, &a[j + j * a_dim1], lda, &work[j], &lbwork, &
00385                     a[j + (*m + 1) * a_dim1], lda, &work[lbwork * nb + nt * 
00386                     nt + 1], &i__3);
00387         } else {
00388             i__2 = *m - j + 1;
00389             i__1 = *n - *m;
00390             i__4 = k - j + 1;
00391             i__3 = *n - *m;
00392             slarfb_("Left", "Transpose", "Forward", "Columnwise", &i__2, &
00393                     i__1, &i__4, &a[j + j * a_dim1], lda, &work[lbwork * nb + 
00394                     1], &nt, &a[j + (*m + 1) * a_dim1], lda, &work[lbwork * 
00395                     nb + nt * nt + 1], &i__3);
00396         }
00397     }
00398     work[1] = (real) iws;
00399     return 0;
00400 
00401 /*     End of SGEQRF */
00402 
00403 } /* sgeqrf_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:06