slaqps.c
Go to the documentation of this file.
00001 /* slaqps.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 real c_b8 = -1.f;
00020 static real c_b9 = 1.f;
00021 static real c_b16 = 0.f;
00022 
00023 /* Subroutine */ int slaqps_(integer *m, integer *n, integer *offset, integer 
00024         *nb, integer *kb, real *a, integer *lda, integer *jpvt, real *tau, 
00025         real *vn1, real *vn2, real *auxv, real *f, integer *ldf)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, f_dim1, f_offset, i__1, i__2;
00029     real r__1, r__2;
00030 
00031     /* Builtin functions */
00032     double sqrt(doublereal);
00033     integer i_nint(real *);
00034 
00035     /* Local variables */
00036     integer j, k, rk;
00037     real akk;
00038     integer pvt;
00039     real temp, temp2;
00040     extern doublereal snrm2_(integer *, real *, integer *);
00041     real tol3z;
00042     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
00043             integer *, real *, real *, integer *, real *, integer *, real *, 
00044             real *, integer *);
00045     integer itemp;
00046     extern /* Subroutine */ int sgemv_(char *, integer *, integer *, real *, 
00047             real *, integer *, real *, integer *, real *, real *, integer *), sswap_(integer *, real *, integer *, real *, integer *);
00048     extern doublereal slamch_(char *);
00049     integer lsticc;
00050     extern integer isamax_(integer *, real *, integer *);
00051     extern /* Subroutine */ int slarfp_(integer *, real *, real *, integer *, 
00052             real *);
00053     integer lastrk;
00054 
00055 
00056 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00057 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00058 /*     November 2006 */
00059 
00060 /*     .. Scalar Arguments .. */
00061 /*     .. */
00062 /*     .. Array Arguments .. */
00063 /*     .. */
00064 
00065 /*  Purpose */
00066 /*  ======= */
00067 
00068 /*  SLAQPS computes a step of QR factorization with column pivoting */
00069 /*  of a real M-by-N matrix A by using Blas-3.  It tries to factorize */
00070 /*  NB columns from A starting from the row OFFSET+1, and updates all */
00071 /*  of the matrix with Blas-3 xGEMM. */
00072 
00073 /*  In some cases, due to catastrophic cancellations, it cannot */
00074 /*  factorize NB columns.  Hence, the actual number of factorized */
00075 /*  columns is returned in KB. */
00076 
00077 /*  Block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. */
00078 
00079 /*  Arguments */
00080 /*  ========= */
00081 
00082 /*  M       (input) INTEGER */
00083 /*          The number of rows of the matrix A. M >= 0. */
00084 
00085 /*  N       (input) INTEGER */
00086 /*          The number of columns of the matrix A. N >= 0 */
00087 
00088 /*  OFFSET  (input) INTEGER */
00089 /*          The number of rows of A that have been factorized in */
00090 /*          previous steps. */
00091 
00092 /*  NB      (input) INTEGER */
00093 /*          The number of columns to factorize. */
00094 
00095 /*  KB      (output) INTEGER */
00096 /*          The number of columns actually factorized. */
00097 
00098 /*  A       (input/output) REAL array, dimension (LDA,N) */
00099 /*          On entry, the M-by-N matrix A. */
00100 /*          On exit, block A(OFFSET+1:M,1:KB) is the triangular */
00101 /*          factor obtained and block A(1:OFFSET,1:N) has been */
00102 /*          accordingly pivoted, but no factorized. */
00103 /*          The rest of the matrix, block A(OFFSET+1:M,KB+1:N) has */
00104 /*          been updated. */
00105 
00106 /*  LDA     (input) INTEGER */
00107 /*          The leading dimension of the array A. LDA >= max(1,M). */
00108 
00109 /*  JPVT    (input/output) INTEGER array, dimension (N) */
00110 /*          JPVT(I) = K <==> Column K of the full matrix A has been */
00111 /*          permuted into position I in AP. */
00112 
00113 /*  TAU     (output) REAL array, dimension (KB) */
00114 /*          The scalar factors of the elementary reflectors. */
00115 
00116 /*  VN1     (input/output) REAL array, dimension (N) */
00117 /*          The vector with the partial column norms. */
00118 
00119 /*  VN2     (input/output) REAL array, dimension (N) */
00120 /*          The vector with the exact column norms. */
00121 
00122 /*  AUXV    (input/output) REAL array, dimension (NB) */
00123 /*          Auxiliar vector. */
00124 
00125 /*  F       (input/output) REAL array, dimension (LDF,NB) */
00126 /*          Matrix F' = L*Y'*A. */
00127 
00128 /*  LDF     (input) INTEGER */
00129 /*          The leading dimension of the array F. LDF >= max(1,N). */
00130 
00131 /*  Further Details */
00132 /*  =============== */
00133 
00134 /*  Based on contributions by */
00135 /*    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain */
00136 /*    X. Sun, Computer Science Dept., Duke University, USA */
00137 
00138 /*  Partial column norm updating strategy modified by */
00139 /*    Z. Drmac and Z. Bujanovic, Dept. of Mathematics, */
00140 /*    University of Zagreb, Croatia. */
00141 /*    June 2006. */
00142 /*  For more details see LAPACK Working Note 176. */
00143 /*  ===================================================================== */
00144 
00145 /*     .. Parameters .. */
00146 /*     .. */
00147 /*     .. Local Scalars .. */
00148 /*     .. */
00149 /*     .. External Subroutines .. */
00150 /*     .. */
00151 /*     .. Intrinsic Functions .. */
00152 /*     .. */
00153 /*     .. External Functions .. */
00154 /*     .. */
00155 /*     .. Executable Statements .. */
00156 
00157     /* Parameter adjustments */
00158     a_dim1 = *lda;
00159     a_offset = 1 + a_dim1;
00160     a -= a_offset;
00161     --jpvt;
00162     --tau;
00163     --vn1;
00164     --vn2;
00165     --auxv;
00166     f_dim1 = *ldf;
00167     f_offset = 1 + f_dim1;
00168     f -= f_offset;
00169 
00170     /* Function Body */
00171 /* Computing MIN */
00172     i__1 = *m, i__2 = *n + *offset;
00173     lastrk = min(i__1,i__2);
00174     lsticc = 0;
00175     k = 0;
00176     tol3z = sqrt(slamch_("Epsilon"));
00177 
00178 /*     Beginning of while loop. */
00179 
00180 L10:
00181     if (k < *nb && lsticc == 0) {
00182         ++k;
00183         rk = *offset + k;
00184 
00185 /*        Determine ith pivot column and swap if necessary */
00186 
00187         i__1 = *n - k + 1;
00188         pvt = k - 1 + isamax_(&i__1, &vn1[k], &c__1);
00189         if (pvt != k) {
00190             sswap_(m, &a[pvt * a_dim1 + 1], &c__1, &a[k * a_dim1 + 1], &c__1);
00191             i__1 = k - 1;
00192             sswap_(&i__1, &f[pvt + f_dim1], ldf, &f[k + f_dim1], ldf);
00193             itemp = jpvt[pvt];
00194             jpvt[pvt] = jpvt[k];
00195             jpvt[k] = itemp;
00196             vn1[pvt] = vn1[k];
00197             vn2[pvt] = vn2[k];
00198         }
00199 
00200 /*        Apply previous Householder reflectors to column K: */
00201 /*        A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'. */
00202 
00203         if (k > 1) {
00204             i__1 = *m - rk + 1;
00205             i__2 = k - 1;
00206             sgemv_("No transpose", &i__1, &i__2, &c_b8, &a[rk + a_dim1], lda, 
00207                     &f[k + f_dim1], ldf, &c_b9, &a[rk + k * a_dim1], &c__1);
00208         }
00209 
00210 /*        Generate elementary reflector H(k). */
00211 
00212         if (rk < *m) {
00213             i__1 = *m - rk + 1;
00214             slarfp_(&i__1, &a[rk + k * a_dim1], &a[rk + 1 + k * a_dim1], &
00215                     c__1, &tau[k]);
00216         } else {
00217             slarfp_(&c__1, &a[rk + k * a_dim1], &a[rk + k * a_dim1], &c__1, &
00218                     tau[k]);
00219         }
00220 
00221         akk = a[rk + k * a_dim1];
00222         a[rk + k * a_dim1] = 1.f;
00223 
00224 /*        Compute Kth column of F: */
00225 
00226 /*        Compute  F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K). */
00227 
00228         if (k < *n) {
00229             i__1 = *m - rk + 1;
00230             i__2 = *n - k;
00231             sgemv_("Transpose", &i__1, &i__2, &tau[k], &a[rk + (k + 1) * 
00232                     a_dim1], lda, &a[rk + k * a_dim1], &c__1, &c_b16, &f[k + 
00233                     1 + k * f_dim1], &c__1);
00234         }
00235 
00236 /*        Padding F(1:K,K) with zeros. */
00237 
00238         i__1 = k;
00239         for (j = 1; j <= i__1; ++j) {
00240             f[j + k * f_dim1] = 0.f;
00241 /* L20: */
00242         }
00243 
00244 /*        Incremental updating of F: */
00245 /*        F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)' */
00246 /*                    *A(RK:M,K). */
00247 
00248         if (k > 1) {
00249             i__1 = *m - rk + 1;
00250             i__2 = k - 1;
00251             r__1 = -tau[k];
00252             sgemv_("Transpose", &i__1, &i__2, &r__1, &a[rk + a_dim1], lda, &a[
00253                     rk + k * a_dim1], &c__1, &c_b16, &auxv[1], &c__1);
00254 
00255             i__1 = k - 1;
00256             sgemv_("No transpose", n, &i__1, &c_b9, &f[f_dim1 + 1], ldf, &
00257                     auxv[1], &c__1, &c_b9, &f[k * f_dim1 + 1], &c__1);
00258         }
00259 
00260 /*        Update the current row of A: */
00261 /*        A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'. */
00262 
00263         if (k < *n) {
00264             i__1 = *n - k;
00265             sgemv_("No transpose", &i__1, &k, &c_b8, &f[k + 1 + f_dim1], ldf, 
00266                     &a[rk + a_dim1], lda, &c_b9, &a[rk + (k + 1) * a_dim1], 
00267                     lda);
00268         }
00269 
00270 /*        Update partial column norms. */
00271 
00272         if (rk < lastrk) {
00273             i__1 = *n;
00274             for (j = k + 1; j <= i__1; ++j) {
00275                 if (vn1[j] != 0.f) {
00276 
00277 /*                 NOTE: The following 4 lines follow from the analysis in */
00278 /*                 Lapack Working Note 176. */
00279 
00280                     temp = (r__1 = a[rk + j * a_dim1], dabs(r__1)) / vn1[j];
00281 /* Computing MAX */
00282                     r__1 = 0.f, r__2 = (temp + 1.f) * (1.f - temp);
00283                     temp = dmax(r__1,r__2);
00284 /* Computing 2nd power */
00285                     r__1 = vn1[j] / vn2[j];
00286                     temp2 = temp * (r__1 * r__1);
00287                     if (temp2 <= tol3z) {
00288                         vn2[j] = (real) lsticc;
00289                         lsticc = j;
00290                     } else {
00291                         vn1[j] *= sqrt(temp);
00292                     }
00293                 }
00294 /* L30: */
00295             }
00296         }
00297 
00298         a[rk + k * a_dim1] = akk;
00299 
00300 /*        End of while loop. */
00301 
00302         goto L10;
00303     }
00304     *kb = k;
00305     rk = *offset + *kb;
00306 
00307 /*     Apply the block reflector to the rest of the matrix: */
00308 /*     A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) - */
00309 /*                         A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'. */
00310 
00311 /* Computing MIN */
00312     i__1 = *n, i__2 = *m - *offset;
00313     if (*kb < min(i__1,i__2)) {
00314         i__1 = *m - rk;
00315         i__2 = *n - *kb;
00316         sgemm_("No transpose", "Transpose", &i__1, &i__2, kb, &c_b8, &a[rk + 
00317                 1 + a_dim1], lda, &f[*kb + 1 + f_dim1], ldf, &c_b9, &a[rk + 1 
00318                 + (*kb + 1) * a_dim1], lda);
00319     }
00320 
00321 /*     Recomputation of difficult columns. */
00322 
00323 L40:
00324     if (lsticc > 0) {
00325         itemp = i_nint(&vn2[lsticc]);
00326         i__1 = *m - rk;
00327         vn1[lsticc] = snrm2_(&i__1, &a[rk + 1 + lsticc * a_dim1], &c__1);
00328 
00329 /*        NOTE: The computation of VN1( LSTICC ) relies on the fact that */
00330 /*        SNRM2 does not fail on vectors with norm below the value of */
00331 /*        SQRT(DLAMCH('S')) */
00332 
00333         vn2[lsticc] = vn1[lsticc];
00334         lsticc = itemp;
00335         goto L40;
00336     }
00337 
00338     return 0;
00339 
00340 /*     End of SLAQPS */
00341 
00342 } /* slaqps_ */


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