swiftnav: zunmbr.c Source File

Go to the documentation of this file.
00001 /* zunmbr.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__2 = 2;
00021 
00022 /* Subroutine */ int zunmbr_(char *vect, char *side, char *trans, integer *m, 
00023         integer *n, integer *k, doublecomplex *a, integer *lda, doublecomplex 
00024         *tau, doublecomplex *c__, integer *ldc, doublecomplex *work, integer *
00025         lwork, integer *info)
00026 {
00027     /* System generated locals */
00028     address a__1[2];
00029     integer a_dim1, a_offset, c_dim1, c_offset, i__1, i__2, i__3[2];
00030     char ch__1[2];
00031 
00032     /* Builtin functions */
00033     /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
00034 
00035     /* Local variables */
00036     integer i1, i2, nb, mi, ni, nq, nw;
00037     logical left;
00038     extern logical lsame_(char *, char *);
00039     integer iinfo;
00040     extern /* Subroutine */ int xerbla_(char *, integer *);
00041     extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 
00042             integer *, integer *);
00043     logical notran, applyq;
00044     char transt[1];
00045     integer lwkopt;
00046     logical lquery;
00047     extern /* Subroutine */ int zunmlq_(char *, char *, integer *, integer *, 
00048             integer *, doublecomplex *, integer *, doublecomplex *, 
00049             doublecomplex *, integer *, doublecomplex *, integer *, integer *), zunmqr_(char *, char *, integer *, integer *, 
00050             integer *, doublecomplex *, integer *, doublecomplex *, 
00051             doublecomplex *, integer *, doublecomplex *, integer *, integer *);
00052 
00053 
00054 /*  -- LAPACK routine (version 3.2) -- */
00055 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00056 /*     November 2006 */
00057 
00058 /*     .. Scalar Arguments .. */
00059 /*     .. */
00060 /*     .. Array Arguments .. */
00061 /*     .. */
00062 
00063 /*  Purpose */
00064 /*  ======= */
00065 
00066 /*  If VECT = 'Q', ZUNMBR overwrites the general complex M-by-N matrix C */
00067 /*  with */
00068 /*                  SIDE = 'L'     SIDE = 'R' */
00069 /*  TRANS = 'N':      Q * C          C * Q */
00070 /*  TRANS = 'C':      Q**H * C       C * Q**H */
00071 
00072 /*  If VECT = 'P', ZUNMBR overwrites the general complex M-by-N matrix C */
00073 /*  with */
00074 /*                  SIDE = 'L'     SIDE = 'R' */
00075 /*  TRANS = 'N':      P * C          C * P */
00076 /*  TRANS = 'C':      P**H * C       C * P**H */
00077 
00078 /*  Here Q and P**H are the unitary matrices determined by ZGEBRD when */
00079 /*  reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q */
00080 /*  and P**H are defined as products of elementary reflectors H(i) and */
00081 /*  G(i) respectively. */
00082 
00083 /*  Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the */
00084 /*  order of the unitary matrix Q or P**H that is applied. */
00085 
00086 /*  If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: */
00087 /*  if nq >= k, Q = H(1) H(2) . . . H(k); */
00088 /*  if nq < k, Q = H(1) H(2) . . . H(nq-1). */
00089 
00090 /*  If VECT = 'P', A is assumed to have been a K-by-NQ matrix: */
00091 /*  if k < nq, P = G(1) G(2) . . . G(k); */
00092 /*  if k >= nq, P = G(1) G(2) . . . G(nq-1). */
00093 
00094 /*  Arguments */
00095 /*  ========= */
00096 
00097 /*  VECT    (input) CHARACTER*1 */
00098 /*          = 'Q': apply Q or Q**H; */
00099 /*          = 'P': apply P or P**H. */
00100 
00101 /*  SIDE    (input) CHARACTER*1 */
00102 /*          = 'L': apply Q, Q**H, P or P**H from the Left; */
00103 /*          = 'R': apply Q, Q**H, P or P**H from the Right. */
00104 
00105 /*  TRANS   (input) CHARACTER*1 */
00106 /*          = 'N':  No transpose, apply Q or P; */
00107 /*          = 'C':  Conjugate transpose, apply Q**H or P**H. */
00108 
00109 /*  M       (input) INTEGER */
00110 /*          The number of rows of the matrix C. M >= 0. */
00111 
00112 /*  N       (input) INTEGER */
00113 /*          The number of columns of the matrix C. N >= 0. */
00114 
00115 /*  K       (input) INTEGER */
00116 /*          If VECT = 'Q', the number of columns in the original */
00117 /*          matrix reduced by ZGEBRD. */
00118 /*          If VECT = 'P', the number of rows in the original */
00119 /*          matrix reduced by ZGEBRD. */
00120 /*          K >= 0. */
00121 
00122 /*  A       (input) COMPLEX*16 array, dimension */
00123 /*                                (LDA,min(nq,K)) if VECT = 'Q' */
00124 /*                                (LDA,nq)        if VECT = 'P' */
00125 /*          The vectors which define the elementary reflectors H(i) and */
00126 /*          G(i), whose products determine the matrices Q and P, as */
00127 /*          returned by ZGEBRD. */
00128 
00129 /*  LDA     (input) INTEGER */
00130 /*          The leading dimension of the array A. */
00131 /*          If VECT = 'Q', LDA >= max(1,nq); */
00132 /*          if VECT = 'P', LDA >= max(1,min(nq,K)). */
00133 
00134 /*  TAU     (input) COMPLEX*16 array, dimension (min(nq,K)) */
00135 /*          TAU(i) must contain the scalar factor of the elementary */
00136 /*          reflector H(i) or G(i) which determines Q or P, as returned */
00137 /*          by ZGEBRD in the array argument TAUQ or TAUP. */
00138 
00139 /*  C       (input/output) COMPLEX*16 array, dimension (LDC,N) */
00140 /*          On entry, the M-by-N matrix C. */
00141 /*          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q */
00142 /*          or P*C or P**H*C or C*P or C*P**H. */
00143 
00144 /*  LDC     (input) INTEGER */
00145 /*          The leading dimension of the array C. LDC >= max(1,M). */
00146 
00147 /*  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) */
00148 /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
00149 
00150 /*  LWORK   (input) INTEGER */
00151 /*          The dimension of the array WORK. */
00152 /*          If SIDE = 'L', LWORK >= max(1,N); */
00153 /*          if SIDE = 'R', LWORK >= max(1,M); */
00154 /*          if N = 0 or M = 0, LWORK >= 1. */
00155 /*          For optimum performance LWORK >= max(1,N*NB) if SIDE = 'L', */
00156 /*          and LWORK >= max(1,M*NB) if SIDE = 'R', where NB is the */
00157 /*          optimal blocksize. (NB = 0 if M = 0 or N = 0.) */
00158 
00159 /*          If LWORK = -1, then a workspace query is assumed; the routine */
00160 /*          only calculates the optimal size of the WORK array, returns */
00161 /*          this value as the first entry of the WORK array, and no error */
00162 /*          message related to LWORK is issued by XERBLA. */
00163 
00164 /*  INFO    (output) INTEGER */
00165 /*          = 0:  successful exit */
00166 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
00167 
00168 /*  ===================================================================== */
00169 
00170 /*     .. Local Scalars .. */
00171 /*     .. */
00172 /*     .. External Functions .. */
00173 /*     .. */
00174 /*     .. External Subroutines .. */
00175 /*     .. */
00176 /*     .. Intrinsic Functions .. */
00177 /*     .. */
00178 /*     .. Executable Statements .. */
00179 
00180 /*     Test the input arguments */
00181 
00182     /* Parameter adjustments */
00183     a_dim1 = *lda;
00184     a_offset = 1 + a_dim1;
00185     a -= a_offset;
00186     --tau;
00187     c_dim1 = *ldc;
00188     c_offset = 1 + c_dim1;
00189     c__ -= c_offset;
00190     --work;
00191 
00192     /* Function Body */
00193     *info = 0;
00194     applyq = lsame_(vect, "Q");
00195     left = lsame_(side, "L");
00196     notran = lsame_(trans, "N");
00197     lquery = *lwork == -1;
00198 
00199 /*     NQ is the order of Q or P and NW is the minimum dimension of WORK */
00200 
00201     if (left) {
00202         nq = *m;
00203         nw = *n;
00204     } else {
00205         nq = *n;
00206         nw = *m;
00207     }
00208     if (*m == 0 || *n == 0) {
00209         nw = 0;
00210     }
00211     if (! applyq && ! lsame_(vect, "P")) {
00212         *info = -1;
00213     } else if (! left && ! lsame_(side, "R")) {
00214         *info = -2;
00215     } else if (! notran && ! lsame_(trans, "C")) {
00216         *info = -3;
00217     } else if (*m < 0) {
00218         *info = -4;
00219     } else if (*n < 0) {
00220         *info = -5;
00221     } else if (*k < 0) {
00222         *info = -6;
00223     } else /* if(complicated condition) */ {
00224 /* Computing MAX */
00225         i__1 = 1, i__2 = min(nq,*k);
00226         if (applyq && *lda < max(1,nq) || ! applyq && *lda < max(i__1,i__2)) {
00227             *info = -8;
00228         } else if (*ldc < max(1,*m)) {
00229             *info = -11;
00230         } else if (*lwork < max(1,nw) && ! lquery) {
00231             *info = -13;
00232         }
00233     }
00234 
00235     if (*info == 0) {
00236         if (nw > 0) {
00237             if (applyq) {
00238                 if (left) {
00239 /* Writing concatenation */
00240                     i__3[0] = 1, a__1[0] = side;
00241                     i__3[1] = 1, a__1[1] = trans;
00242                     s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
00243                     i__1 = *m - 1;
00244                     i__2 = *m - 1;
00245                     nb = ilaenv_(&c__1, "ZUNMQR", ch__1, &i__1, n, &i__2, &
00246                             c_n1);
00247                 } else {
00248 /* Writing concatenation */
00249                     i__3[0] = 1, a__1[0] = side;
00250                     i__3[1] = 1, a__1[1] = trans;
00251                     s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
00252                     i__1 = *n - 1;
00253                     i__2 = *n - 1;
00254                     nb = ilaenv_(&c__1, "ZUNMQR", ch__1, m, &i__1, &i__2, &
00255                             c_n1);
00256                 }
00257             } else {
00258                 if (left) {
00259 /* Writing concatenation */
00260                     i__3[0] = 1, a__1[0] = side;
00261                     i__3[1] = 1, a__1[1] = trans;
00262                     s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
00263                     i__1 = *m - 1;
00264                     i__2 = *m - 1;
00265                     nb = ilaenv_(&c__1, "ZUNMLQ", ch__1, &i__1, n, &i__2, &
00266                             c_n1);
00267                 } else {
00268 /* Writing concatenation */
00269                     i__3[0] = 1, a__1[0] = side;
00270                     i__3[1] = 1, a__1[1] = trans;
00271                     s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2);
00272                     i__1 = *n - 1;
00273                     i__2 = *n - 1;
00274                     nb = ilaenv_(&c__1, "ZUNMLQ", ch__1, m, &i__1, &i__2, &
00275                             c_n1);
00276                 }
00277             }
00278 /* Computing MAX */
00279             i__1 = 1, i__2 = nw * nb;
00280             lwkopt = max(i__1,i__2);
00281         } else {
00282             lwkopt = 1;
00283         }
00284         work[1].r = (doublereal) lwkopt, work[1].i = 0.;
00285     }
00286 
00287     if (*info != 0) {
00288         i__1 = -(*info);
00289         xerbla_("ZUNMBR", &i__1);
00290         return 0;
00291     } else if (lquery) {
00292         return 0;
00293     }
00294 
00295 /*     Quick return if possible */
00296 
00297     if (*m == 0 || *n == 0) {
00298         return 0;
00299     }
00300 
00301     if (applyq) {
00302 
00303 /*        Apply Q */
00304 
00305         if (nq >= *k) {
00306 
00307 /*           Q was determined by a call to ZGEBRD with nq >= k */
00308 
00309             zunmqr_(side, trans, m, n, k, &a[a_offset], lda, &tau[1], &c__[
00310                     c_offset], ldc, &work[1], lwork, &iinfo);
00311         } else if (nq > 1) {
00312 
00313 /*           Q was determined by a call to ZGEBRD with nq < k */
00314 
00315             if (left) {
00316                 mi = *m - 1;
00317                 ni = *n;
00318                 i1 = 2;
00319                 i2 = 1;
00320             } else {
00321                 mi = *m;
00322                 ni = *n - 1;
00323                 i1 = 1;
00324                 i2 = 2;
00325             }
00326             i__1 = nq - 1;
00327             zunmqr_(side, trans, &mi, &ni, &i__1, &a[a_dim1 + 2], lda, &tau[1]
00328 , &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &iinfo);
00329         }
00330     } else {
00331 
00332 /*        Apply P */
00333 
00334         if (notran) {
00335             *(unsigned char *)transt = 'C';
00336         } else {
00337             *(unsigned char *)transt = 'N';
00338         }
00339         if (nq > *k) {
00340 
00341 /*           P was determined by a call to ZGEBRD with nq > k */
00342 
00343             zunmlq_(side, transt, m, n, k, &a[a_offset], lda, &tau[1], &c__[
00344                     c_offset], ldc, &work[1], lwork, &iinfo);
00345         } else if (nq > 1) {
00346 
00347 /*           P was determined by a call to ZGEBRD with nq <= k */
00348 
00349             if (left) {
00350                 mi = *m - 1;
00351                 ni = *n;
00352                 i1 = 2;
00353                 i2 = 1;
00354             } else {
00355                 mi = *m;
00356                 ni = *n - 1;
00357                 i1 = 1;
00358                 i2 = 2;
00359             }
00360             i__1 = nq - 1;
00361             zunmlq_(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda, 
00362                      &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, &
00363                     iinfo);
00364         }
00365     }
00366     work[1].r = (doublereal) lwkopt, work[1].i = 0.;
00367     return 0;
00368 
00369 /*     End of ZUNMBR */
00370 
00371 } /* zunmbr_ */