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


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