dggqrf.c
Go to the documentation of this file.
00001 /* dggqrf.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 
00021 /* Subroutine */ int dggqrf_(integer *n, integer *m, integer *p, doublereal *
00022         a, integer *lda, doublereal *taua, doublereal *b, integer *ldb, 
00023         doublereal *taub, doublereal *work, integer *lwork, integer *info)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2;
00027 
00028     /* Local variables */
00029     integer nb, nb1, nb2, nb3, lopt;
00030     extern /* Subroutine */ int dgeqrf_(integer *, integer *, doublereal *, 
00031             integer *, doublereal *, doublereal *, integer *, integer *), 
00032             dgerqf_(integer *, integer *, doublereal *, integer *, doublereal 
00033             *, doublereal *, integer *, integer *), xerbla_(char *, integer *);
00034     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00035             integer *, integer *);
00036     extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, 
00037             integer *, doublereal *, integer *, doublereal *, doublereal *, 
00038             integer *, doublereal *, integer *, integer *);
00039     integer lwkopt;
00040     logical lquery;
00041 
00042 
00043 /*  -- LAPACK routine (version 3.2) -- */
00044 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00045 /*     November 2006 */
00046 
00047 /*     .. Scalar Arguments .. */
00048 /*     .. */
00049 /*     .. Array Arguments .. */
00050 /*     .. */
00051 
00052 /*  Purpose */
00053 /*  ======= */
00054 
00055 /*  DGGQRF computes a generalized QR factorization of an N-by-M matrix A */
00056 /*  and an N-by-P matrix B: */
00057 
00058 /*              A = Q*R,        B = Q*T*Z, */
00059 
00060 /*  where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal */
00061 /*  matrix, and R and T assume one of the forms: */
00062 
00063 /*  if N >= M,  R = ( R11 ) M  ,   or if N < M,  R = ( R11  R12 ) N, */
00064 /*                  (  0  ) N-M                         N   M-N */
00065 /*                     M */
00066 
00067 /*  where R11 is upper triangular, and */
00068 
00069 /*  if N <= P,  T = ( 0  T12 ) N,   or if N > P,  T = ( T11 ) N-P, */
00070 /*                   P-N  N                           ( T21 ) P */
00071 /*                                                       P */
00072 
00073 /*  where T12 or T21 is upper triangular. */
00074 
00075 /*  In particular, if B is square and nonsingular, the GQR factorization */
00076 /*  of A and B implicitly gives the QR factorization of inv(B)*A: */
00077 
00078 /*               inv(B)*A = Z'*(inv(T)*R) */
00079 
00080 /*  where inv(B) denotes the inverse of the matrix B, and Z' denotes the */
00081 /*  transpose of the matrix Z. */
00082 
00083 /*  Arguments */
00084 /*  ========= */
00085 
00086 /*  N       (input) INTEGER */
00087 /*          The number of rows of the matrices A and B. N >= 0. */
00088 
00089 /*  M       (input) INTEGER */
00090 /*          The number of columns of the matrix A.  M >= 0. */
00091 
00092 /*  P       (input) INTEGER */
00093 /*          The number of columns of the matrix B.  P >= 0. */
00094 
00095 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,M) */
00096 /*          On entry, the N-by-M matrix A. */
00097 /*          On exit, the elements on and above the diagonal of the array */
00098 /*          contain the min(N,M)-by-M upper trapezoidal matrix R (R is */
00099 /*          upper triangular if N >= M); the elements below the diagonal, */
00100 /*          with the array TAUA, represent the orthogonal matrix Q as a */
00101 /*          product of min(N,M) elementary reflectors (see Further */
00102 /*          Details). */
00103 
00104 /*  LDA     (input) INTEGER */
00105 /*          The leading dimension of the array A. LDA >= max(1,N). */
00106 
00107 /*  TAUA    (output) DOUBLE PRECISION array, dimension (min(N,M)) */
00108 /*          The scalar factors of the elementary reflectors which */
00109 /*          represent the orthogonal matrix Q (see Further Details). */
00110 
00111 /*  B       (input/output) DOUBLE PRECISION array, dimension (LDB,P) */
00112 /*          On entry, the N-by-P matrix B. */
00113 /*          On exit, if N <= P, the upper triangle of the subarray */
00114 /*          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; */
00115 /*          if N > P, the elements on and above the (N-P)-th subdiagonal */
00116 /*          contain the N-by-P upper trapezoidal matrix T; the remaining */
00117 /*          elements, with the array TAUB, represent the orthogonal */
00118 /*          matrix Z as a product of elementary reflectors (see Further */
00119 /*          Details). */
00120 
00121 /*  LDB     (input) INTEGER */
00122 /*          The leading dimension of the array B. LDB >= max(1,N). */
00123 
00124 /*  TAUB    (output) DOUBLE PRECISION array, dimension (min(N,P)) */
00125 /*          The scalar factors of the elementary reflectors which */
00126 /*          represent the orthogonal matrix Z (see Further Details). */
00127 
00128 /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
00129 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00130 
00131 /*  LWORK   (input) INTEGER */
00132 /*          The dimension of the array WORK. LWORK >= max(1,N,M,P). */
00133 /*          For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3), */
00134 /*          where NB1 is the optimal blocksize for the QR factorization */
00135 /*          of an N-by-M matrix, NB2 is the optimal blocksize for the */
00136 /*          RQ factorization of an N-by-P matrix, and NB3 is the optimal */
00137 /*          blocksize for a call of DORMQR. */
00138 
00139 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00140 /*          only calculates the optimal size of the WORK array, returns */
00141 /*          this value as the first entry of the WORK array, and no error */
00142 /*          message related to LWORK is issued by XERBLA. */
00143 
00144 /*  INFO    (output) INTEGER */
00145 /*          = 0:  successful exit */
00146 /*          < 0:  if INFO = -i, the i-th argument had an illegal value. */
00147 
00148 /*  Further Details */
00149 /*  =============== */
00150 
00151 /*  The matrix Q is represented as a product of elementary reflectors */
00152 
00153 /*     Q = H(1) H(2) . . . H(k), where k = min(n,m). */
00154 
00155 /*  Each H(i) has the form */
00156 
00157 /*     H(i) = I - taua * v * v' */
00158 
00159 /*  where taua is a real scalar, and v is a real vector with */
00160 /*  v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), */
00161 /*  and taua in TAUA(i). */
00162 /*  To form Q explicitly, use LAPACK subroutine DORGQR. */
00163 /*  To use Q to update another matrix, use LAPACK subroutine DORMQR. */
00164 
00165 /*  The matrix Z is represented as a product of elementary reflectors */
00166 
00167 /*     Z = H(1) H(2) . . . H(k), where k = min(n,p). */
00168 
00169 /*  Each H(i) has the form */
00170 
00171 /*     H(i) = I - taub * v * v' */
00172 
00173 /*  where taub is a real scalar, and v is a real vector with */
00174 /*  v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in */
00175 /*  B(n-k+i,1:p-k+i-1), and taub in TAUB(i). */
00176 /*  To form Z explicitly, use LAPACK subroutine DORGRQ. */
00177 /*  To use Z to update another matrix, use LAPACK subroutine DORMRQ. */
00178 
00179 /*  ===================================================================== */
00180 
00181 /*     .. Local Scalars .. */
00182 /*     .. */
00183 /*     .. External Subroutines .. */
00184 /*     .. */
00185 /*     .. External Functions .. */
00186 /*     .. */
00187 /*     .. Intrinsic Functions .. */
00188 /*     .. */
00189 /*     .. Executable Statements .. */
00190 
00191 /*     Test the input parameters */
00192 
00193     /* Parameter adjustments */
00194     a_dim1 = *lda;
00195     a_offset = 1 + a_dim1;
00196     a -= a_offset;
00197     --taua;
00198     b_dim1 = *ldb;
00199     b_offset = 1 + b_dim1;
00200     b -= b_offset;
00201     --taub;
00202     --work;
00203 
00204     /* Function Body */
00205     *info = 0;
00206     nb1 = ilaenv_(&c__1, "DGEQRF", " ", n, m, &c_n1, &c_n1);
00207     nb2 = ilaenv_(&c__1, "DGERQF", " ", n, p, &c_n1, &c_n1);
00208     nb3 = ilaenv_(&c__1, "DORMQR", " ", n, m, p, &c_n1);
00209 /* Computing MAX */
00210     i__1 = max(nb1,nb2);
00211     nb = max(i__1,nb3);
00212 /* Computing MAX */
00213     i__1 = max(*n,*m);
00214     lwkopt = max(i__1,*p) * nb;
00215     work[1] = (doublereal) lwkopt;
00216     lquery = *lwork == -1;
00217     if (*n < 0) {
00218         *info = -1;
00219     } else if (*m < 0) {
00220         *info = -2;
00221     } else if (*p < 0) {
00222         *info = -3;
00223     } else if (*lda < max(1,*n)) {
00224         *info = -5;
00225     } else if (*ldb < max(1,*n)) {
00226         *info = -8;
00227     } else /* if(complicated condition) */ {
00228 /* Computing MAX */
00229         i__1 = max(1,*n), i__1 = max(i__1,*m);
00230         if (*lwork < max(i__1,*p) && ! lquery) {
00231             *info = -11;
00232         }
00233     }
00234     if (*info != 0) {
00235         i__1 = -(*info);
00236         xerbla_("DGGQRF", &i__1);
00237         return 0;
00238     } else if (lquery) {
00239         return 0;
00240     }
00241 
00242 /*     QR factorization of N-by-M matrix A: A = Q*R */
00243 
00244     dgeqrf_(n, m, &a[a_offset], lda, &taua[1], &work[1], lwork, info);
00245     lopt = (integer) work[1];
00246 
00247 /*     Update B := Q'*B. */
00248 
00249     i__1 = min(*n,*m);
00250     dormqr_("Left", "Transpose", n, p, &i__1, &a[a_offset], lda, &taua[1], &b[
00251             b_offset], ldb, &work[1], lwork, info);
00252 /* Computing MAX */
00253     i__1 = lopt, i__2 = (integer) work[1];
00254     lopt = max(i__1,i__2);
00255 
00256 /*     RQ factorization of N-by-P matrix B: B = T*Z. */
00257 
00258     dgerqf_(n, p, &b[b_offset], ldb, &taub[1], &work[1], lwork, info);
00259 /* Computing MAX */
00260     i__1 = lopt, i__2 = (integer) work[1];
00261     work[1] = (doublereal) max(i__1,i__2);
00262 
00263     return 0;
00264 
00265 /*     End of DGGQRF */
00266 
00267 } /* dggqrf_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:45