/* sggrqf.f -- translated by f2c (version 20061008).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
#include "blaswrap.h"
/* Table of constant values */
static integer c__1 = 1;
static integer c_n1 = -1;
/* Subroutine */ int sggrqf_(integer *m, integer *p, integer *n, real *a,
integer *lda, real *taua, real *b, integer *ldb, real *taub, real *
work, integer *lwork, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3;
/* Local variables */
integer nb, nb1, nb2, nb3, lopt;
extern /* Subroutine */ int xerbla_(char *, integer *);
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
integer *, integer *);
extern /* Subroutine */ int sgeqrf_(integer *, integer *, real *, integer
*, real *, real *, integer *, integer *), sgerqf_(integer *,
integer *, real *, integer *, real *, real *, integer *, integer *
);
integer lwkopt;
logical lquery;
extern /* Subroutine */ int sormrq_(char *, char *, integer *, integer *,
integer *, real *, integer *, real *, real *, integer *, real *,
integer *, integer *);
/* -- LAPACK routine (version 3.2) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* SGGRQF computes a generalized RQ factorization of an M-by-N matrix A */
/* and a P-by-N matrix B: */
/* A = R*Q, B = Z*T*Q, */
/* where Q is an N-by-N orthogonal matrix, Z is a P-by-P orthogonal */
/* matrix, and R and T assume one of the forms: */
/* if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N, */
/* N-M M ( R21 ) N */
/* N */
/* where R12 or R21 is upper triangular, and */
/* if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P, */
/* ( 0 ) P-N P N-P */
/* N */
/* where T11 is upper triangular. */
/* In particular, if B is square and nonsingular, the GRQ factorization */
/* of A and B implicitly gives the RQ factorization of A*inv(B): */
/* A*inv(B) = (R*inv(T))*Z' */
/* where inv(B) denotes the inverse of the matrix B, and Z' denotes the */
/* transpose of the matrix Z. */
/* Arguments */
/* ========= */
/* M (input) INTEGER */
/* The number of rows of the matrix A. M >= 0. */
/* P (input) INTEGER */
/* The number of rows of the matrix B. P >= 0. */
/* N (input) INTEGER */
/* The number of columns of the matrices A and B. N >= 0. */
/* A (input/output) REAL array, dimension (LDA,N) */
/* On entry, the M-by-N matrix A. */
/* On exit, if M <= N, the upper triangle of the subarray */
/* A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R; */
/* if M > N, the elements on and above the (M-N)-th subdiagonal */
/* contain the M-by-N upper trapezoidal matrix R; the remaining */
/* elements, with the array TAUA, represent the orthogonal */
/* matrix Q as a product of elementary reflectors (see Further */
/* Details). */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,M). */
/* TAUA (output) REAL array, dimension (min(M,N)) */
/* The scalar factors of the elementary reflectors which */
/* represent the orthogonal matrix Q (see Further Details). */
/* B (input/output) REAL array, dimension (LDB,N) */
/* On entry, the P-by-N matrix B. */
/* On exit, the elements on and above the diagonal of the array */
/* contain the min(P,N)-by-N upper trapezoidal matrix T (T is */
/* upper triangular if P >= N); the elements below the diagonal, */
/* with the array TAUB, represent the orthogonal matrix Z as a */
/* product of elementary reflectors (see Further Details). */
/* LDB (input) INTEGER */
/* The leading dimension of the array B. LDB >= max(1,P). */
/* TAUB (output) REAL array, dimension (min(P,N)) */
/* The scalar factors of the elementary reflectors which */
/* represent the orthogonal matrix Z (see Further Details). */
/* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */
/* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
/* LWORK (input) INTEGER */
/* The dimension of the array WORK. LWORK >= max(1,N,M,P). */
/* For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3), */
/* where NB1 is the optimal blocksize for the RQ factorization */
/* of an M-by-N matrix, NB2 is the optimal blocksize for the */
/* QR factorization of a P-by-N matrix, and NB3 is the optimal */
/* blocksize for a call of SORMRQ. */
/* If LWORK = -1, then a workspace query is assumed; the routine */
/* only calculates the optimal size of the WORK array, returns */
/* this value as the first entry of the WORK array, and no error */
/* message related to LWORK is issued by XERBLA. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INF0= -i, the i-th argument had an illegal value. */
/* Further Details */
/* =============== */
/* The matrix Q is represented as a product of elementary reflectors */
/* Q = H(1) H(2) . . . H(k), where k = min(m,n). */
/* Each H(i) has the form */
/* H(i) = I - taua * v * v' */
/* where taua is a real scalar, and v is a real vector with */
/* v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in */
/* A(m-k+i,1:n-k+i-1), and taua in TAUA(i). */
/* To form Q explicitly, use LAPACK subroutine SORGRQ. */
/* To use Q to update another matrix, use LAPACK subroutine SORMRQ. */
/* The matrix Z is represented as a product of elementary reflectors */
/* Z = H(1) H(2) . . . H(k), where k = min(p,n). */
/* Each H(i) has the form */
/* H(i) = I - taub * v * v' */
/* where taub is a real scalar, and v is a real vector with */
/* v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i), */
/* and taub in TAUB(i). */
/* To form Z explicitly, use LAPACK subroutine SORGQR. */
/* To use Z to update another matrix, use LAPACK subroutine SORMQR. */
/* ===================================================================== */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--taua;
b_dim1 = *ldb;
b_offset = 1 + b_dim1;
b -= b_offset;
--taub;
--work;
/* Function Body */
*info = 0;
nb1 = ilaenv_(&c__1, "SGERQF", " ", m, n, &c_n1, &c_n1);
nb2 = ilaenv_(&c__1, "SGEQRF", " ", p, n, &c_n1, &c_n1);
nb3 = ilaenv_(&c__1, "SORMRQ", " ", m, n, p, &c_n1);
/* Computing MAX */
i__1 = max(nb1,nb2);
nb = max(i__1,nb3);
/* Computing MAX */
i__1 = max(*n,*m);
lwkopt = max(i__1,*p) * nb;
work[1] = (real) lwkopt;
lquery = *lwork == -1;
if (*m < 0) {
*info = -1;
} else if (*p < 0) {
*info = -2;
} else if (*n < 0) {
*info = -3;
} else if (*lda < max(1,*m)) {
*info = -5;
} else if (*ldb < max(1,*p)) {
*info = -8;
} else /* if(complicated condition) */ {
/* Computing MAX */
i__1 = max(1,*m), i__1 = max(i__1,*p);
if (*lwork < max(i__1,*n) && ! lquery) {
*info = -11;
}
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("SGGRQF", &i__1);
return 0;
} else if (lquery) {
return 0;
}
/* RQ factorization of M-by-N matrix A: A = R*Q */
sgerqf_(m, n, &a[a_offset], lda, &taua[1], &work[1], lwork, info);
lopt = work[1];
/* Update B := B*Q' */
i__1 = min(*m,*n);
/* Computing MAX */
i__2 = 1, i__3 = *m - *n + 1;
sormrq_("Right", "Transpose", p, n, &i__1, &a[max(i__2, i__3)+ a_dim1],
lda, &taua[1], &b[b_offset], ldb, &work[1], lwork, info);
/* Computing MAX */
i__1 = lopt, i__2 = (integer) work[1];
lopt = max(i__1,i__2);
/* QR factorization of P-by-N matrix B: B = Z*T */
sgeqrf_(p, n, &b[b_offset], ldb, &taub[1], &work[1], lwork, info);
/* Computing MAX */
i__1 = lopt, i__2 = (integer) work[1];
work[1] = (real) max(i__1,i__2);
return 0;
/* End of SGGRQF */
} /* sggrqf_ */