swiftnav: cunmbr.c Source File

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