cdrvbd.c
Go to the documentation of this file.
00001 /* cdrvbd.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 complex c_b1 = {0.f,0.f};
00019 static complex c_b2 = {1.f,0.f};
00020 static integer c__4 = 4;
00021 static integer c__1 = 1;
00022 static integer c__0 = 0;
00023 
00024 /* Subroutine */ int cdrvbd_(integer *nsizes, integer *mm, integer *nn, 
00025         integer *ntypes, logical *dotype, integer *iseed, real *thresh, 
00026         complex *a, integer *lda, complex *u, integer *ldu, complex *vt, 
00027         integer *ldvt, complex *asav, complex *usav, complex *vtsav, real *s, 
00028         real *ssav, real *e, complex *work, integer *lwork, real *rwork, 
00029         integer *iwork, integer *nounit, integer *info)
00030 {
00031     /* Initialized data */
00032 
00033     static char cjob[1*4] = "N" "O" "S" "A";
00034 
00035     /* Format strings */
00036     static char fmt_9996[] = "(\002 CDRVBD: \002,a,\002 returned INFO=\002,i"
00037             "6,\002.\002,/9x,\002M=\002,i6,\002, N=\002,i6,\002, JTYPE=\002,i"
00038             "6,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00039     static char fmt_9995[] = "(\002 CDRVBD: \002,a,\002 returned INFO=\002,i"
00040             "6,\002.\002,/9x,\002M=\002,i6,\002, N=\002,i6,\002, JTYPE=\002,i"
00041             "6,\002, LSWORK=\002,i6,/9x,\002ISEED=(\002,3(i5,\002,\002),i5"
00042             ",\002)\002)";
00043     static char fmt_9999[] = "(\002 SVD -- Complex Singular Value Decomposit"
00044             "ion Driver \002,/\002 Matrix types (see CDRVBD for details):\002"
00045             ",//\002 1 = Zero matrix\002,/\002 2 = Identity matrix\002,/\002 "
00046             "3 = Evenly spaced singular values near 1\002,/\002 4 = Evenly sp"
00047             "aced singular values near underflow\002,/\002 5 = Evenly spaced "
00048             "singular values near overflow\002,//\002 Tests performed: ( A is"
00049             " dense, U and V are unitary,\002,/19x,\002 S is an array, and Up"
00050             "artial, VTpartial, and\002,/19x,\002 Spartial are partially comp"
00051             "uted U, VT and S),\002,/)";
00052     static char fmt_9998[] = "(\002 Tests performed with Test Threshold ="
00053             " \002,f8.2,/\002 CGESVD: \002,/\002 1 = | A - U diag(S) VT | / ("
00054             " |A| max(M,N) ulp ) \002,/\002 2 = | I - U**T U | / ( M ulp )"
00055             " \002,/\002 3 = | I - VT VT**T | / ( N ulp ) \002,/\002 4 = 0 if"
00056             " S contains min(M,N) nonnegative values in\002,\002 decreasing o"
00057             "rder, else 1/ulp\002,/\002 5 = | U - Upartial | / ( M ulp )\002,/"
00058             "\002 6 = | VT - VTpartial | / ( N ulp )\002,/\002 7 = | S - Spar"
00059             "tial | / ( min(M,N) ulp |S| )\002,/\002 CGESDD: \002,/\002 8 = |"
00060             " A - U diag(S) VT | / ( |A| max(M,N) ulp ) \002,/\002 9 = | I - "
00061             "U**T U | / ( M ulp ) \002,/\00210 = | I - VT VT**T | / ( N ulp ) "
00062             "\002,/\00211 = 0 if S contains min(M,N) nonnegative values in"
00063             "\002,\002 decreasing order, else 1/ulp\002,/\00212 = | U - Upart"
00064             "ial | / ( M ulp )\002,/\00213 = | VT - VTpartial | / ( N ulp "
00065             ")\002,/\00214 = | S - Spartial | / ( min(M,N) ulp |S| )\002,//)";
00066     static char fmt_9997[] = "(\002 M=\002,i5,\002, N=\002,i5,\002, type "
00067             "\002,i1,\002, IWS=\002,i1,\002, seed=\002,4(i4,\002,\002),\002 t"
00068             "est(\002,i1,\002)=\002,g11.4)";
00069 
00070     /* System generated locals */
00071     integer a_dim1, a_offset, asav_dim1, asav_offset, u_dim1, u_offset, 
00072             usav_dim1, usav_offset, vt_dim1, vt_offset, vtsav_dim1, 
00073             vtsav_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, 
00074             i__9, i__10, i__11, i__12, i__13, i__14;
00075     real r__1, r__2, r__3;
00076 
00077     /* Builtin functions */
00078     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00079 
00080     /* Local variables */
00081     integer i__, j, m, n;
00082     real dif, div;
00083     integer ijq, iju;
00084     real ulp;
00085     char jobq[1], jobu[1];
00086     integer mmax, nmax;
00087     real unfl, ovfl;
00088     integer ijvt;
00089     extern /* Subroutine */ int cbdt01_(integer *, integer *, integer *, 
00090             complex *, integer *, complex *, integer *, real *, real *, 
00091             complex *, integer *, complex *, real *, real *);
00092     logical badmm, badnn;
00093     integer nfail, iinfo;
00094     extern /* Subroutine */ int cunt01_(char *, integer *, integer *, complex 
00095             *, integer *, complex *, integer *, real *, real *);
00096     real anorm;
00097     extern /* Subroutine */ int cunt03_(char *, integer *, integer *, integer 
00098             *, integer *, complex *, integer *, complex *, integer *, complex 
00099             *, integer *, real *, real *, integer *);
00100     integer mnmin, mnmax;
00101     char jobvt[1];
00102     integer iwspc, jsize, nerrs, jtype, ntest, iwtmp;
00103     extern /* Subroutine */ int cgesdd_(char *, integer *, integer *, complex 
00104             *, integer *, real *, complex *, integer *, complex *, integer *, 
00105             complex *, integer *, real *, integer *, integer *);
00106     extern doublereal slamch_(char *);
00107     extern /* Subroutine */ int cgesvd_(char *, char *, integer *, integer *, 
00108             complex *, integer *, real *, complex *, integer *, complex *, 
00109             integer *, complex *, integer *, real *, integer *), clacpy_(char *, integer *, integer *, complex *, integer 
00110             *, complex *, integer *), claset_(char *, integer *, 
00111             integer *, complex *, complex *, complex *, integer *);
00112     integer ioldsd[4];
00113     extern /* Subroutine */ int xerbla_(char *, integer *), alasvm_(
00114             char *, integer *, integer *, integer *, integer *), 
00115             clatms_(integer *, integer *, char *, integer *, char *, real *, 
00116             integer *, real *, real *, integer *, integer *, char *, complex *
00117 , integer *, complex *, integer *);
00118     integer ntestf, minwrk;
00119     real ulpinv, result[14];
00120     integer lswork, mtypes, ntestt;
00121 
00122     /* Fortran I/O blocks */
00123     static cilist io___27 = { 0, 0, 0, fmt_9996, 0 };
00124     static cilist io___32 = { 0, 0, 0, fmt_9995, 0 };
00125     static cilist io___39 = { 0, 0, 0, fmt_9995, 0 };
00126     static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };
00127     static cilist io___44 = { 0, 0, 0, fmt_9998, 0 };
00128     static cilist io___45 = { 0, 0, 0, fmt_9997, 0 };
00129 
00130 
00131 
00132 /*  -- LAPACK test routine (version 3.1) -- */
00133 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00134 /*     November 2006 */
00135 
00136 /*     .. Scalar Arguments .. */
00137 /*     .. */
00138 /*     .. Array Arguments .. */
00139 /*     .. */
00140 
00141 /*  Purpose */
00142 /*  ======= */
00143 
00144 /*  CDRVBD checks the singular value decomposition (SVD) driver CGESVD */
00145 /*  and CGESDD. */
00146 /*  CGESVD and CGESDD factors A = U diag(S) VT, where U and VT are */
00147 /*  unitary and diag(S) is diagonal with the entries of the array S on */
00148 /*  its diagonal. The entries of S are the singular values, nonnegative */
00149 /*  and stored in decreasing order.  U and VT can be optionally not */
00150 /*  computed, overwritten on A, or computed partially. */
00151 
00152 /*  A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN. */
00153 /*  U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N. */
00154 
00155 /*  When CDRVBD is called, a number of matrix "sizes" (M's and N's) */
00156 /*  and a number of matrix "types" are specified.  For each size (M,N) */
00157 /*  and each type of matrix, and for the minimal workspace as well as */
00158 /*  workspace adequate to permit blocking, an  M x N  matrix "A" will be */
00159 /*  generated and used to test the SVD routines.  For each matrix, A will */
00160 /*  be factored as A = U diag(S) VT and the following 12 tests computed: */
00161 
00162 /*  Test for CGESVD: */
00163 
00164 /*  (1)   | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */
00165 
00166 /*  (2)   | I - U'U | / ( M ulp ) */
00167 
00168 /*  (3)   | I - VT VT' | / ( N ulp ) */
00169 
00170 /*  (4)   S contains MNMIN nonnegative values in decreasing order. */
00171 /*        (Return 0 if true, 1/ULP if false.) */
00172 
00173 /*  (5)   | U - Upartial | / ( M ulp ) where Upartial is a partially */
00174 /*        computed U. */
00175 
00176 /*  (6)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially */
00177 /*        computed VT. */
00178 
00179 /*  (7)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the */
00180 /*        vector of singular values from the partial SVD */
00181 
00182 /*  Test for CGESDD: */
00183 
00184 /*  (1)   | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */
00185 
00186 /*  (2)   | I - U'U | / ( M ulp ) */
00187 
00188 /*  (3)   | I - VT VT' | / ( N ulp ) */
00189 
00190 /*  (4)   S contains MNMIN nonnegative values in decreasing order. */
00191 /*        (Return 0 if true, 1/ULP if false.) */
00192 
00193 /*  (5)   | U - Upartial | / ( M ulp ) where Upartial is a partially */
00194 /*        computed U. */
00195 
00196 /*  (6)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially */
00197 /*        computed VT. */
00198 
00199 /*  (7)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the */
00200 /*        vector of singular values from the partial SVD */
00201 
00202 /*  The "sizes" are specified by the arrays MM(1:NSIZES) and */
00203 /*  NN(1:NSIZES); the value of each element pair (MM(j),NN(j)) */
00204 /*  specifies one size.  The "types" are specified by a logical array */
00205 /*  DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j" */
00206 /*  will be generated. */
00207 /*  Currently, the list of possible types is: */
00208 
00209 /*  (1)  The zero matrix. */
00210 /*  (2)  The identity matrix. */
00211 /*  (3)  A matrix of the form  U D V, where U and V are unitary and */
00212 /*       D has evenly spaced entries 1, ..., ULP with random signs */
00213 /*       on the diagonal. */
00214 /*  (4)  Same as (3), but multiplied by the underflow-threshold / ULP. */
00215 /*  (5)  Same as (3), but multiplied by the overflow-threshold * ULP. */
00216 
00217 /*  Arguments */
00218 /*  ========== */
00219 
00220 /*  NSIZES  (input) INTEGER */
00221 /*          The number of sizes of matrices to use.  If it is zero, */
00222 /*          CDRVBD does nothing.  It must be at least zero. */
00223 
00224 /*  MM      (input) INTEGER array, dimension (NSIZES) */
00225 /*          An array containing the matrix "heights" to be used.  For */
00226 /*          each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j) */
00227 /*          will be ignored.  The MM(j) values must be at least zero. */
00228 
00229 /*  NN      (input) INTEGER array, dimension (NSIZES) */
00230 /*          An array containing the matrix "widths" to be used.  For */
00231 /*          each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j) */
00232 /*          will be ignored.  The NN(j) values must be at least zero. */
00233 
00234 /*  NTYPES  (input) INTEGER */
00235 /*          The number of elements in DOTYPE.   If it is zero, CDRVBD */
00236 /*          does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00237 /*          and NSIZES is 1, then an additional type, MAXTYP+1 is */
00238 /*          defined, which is to use whatever matrices are in A and B. */
00239 /*          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00240 /*          DOTYPE(MAXTYP+1) is .TRUE. . */
00241 
00242 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00243 /*          If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix */
00244 /*          of type j will be generated.  If NTYPES is smaller than the */
00245 /*          maximum number of types defined (PARAMETER MAXTYP), then */
00246 /*          types NTYPES+1 through MAXTYP will not be generated.  If */
00247 /*          NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through */
00248 /*          DOTYPE(NTYPES) will be ignored. */
00249 
00250 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00251 /*          On entry ISEED specifies the seed of the random number */
00252 /*          generator. The array elements should be between 0 and 4095; */
00253 /*          if not they will be reduced mod 4096.  Also, ISEED(4) must */
00254 /*          be odd.  The random number generator uses a linear */
00255 /*          congruential sequence limited to small integers, and so */
00256 /*          should produce machine independent random numbers. The */
00257 /*          values of ISEED are changed on exit, and can be used in the */
00258 /*          next call to CDRVBD to continue the same random number */
00259 /*          sequence. */
00260 
00261 /*  THRESH  (input) REAL */
00262 /*          A test will count as "failed" if the "error", computed as */
00263 /*          described above, exceeds THRESH.  Note that the error */
00264 /*          is scaled to be O(1), so THRESH should be a reasonably */
00265 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00266 /*          it should not depend on the precision (single vs. double) */
00267 /*          or the size of the matrix.  It must be at least zero. */
00268 
00269 /*  NOUNIT  (input) INTEGER */
00270 /*          The FORTRAN unit number for printing out error messages */
00271 /*          (e.g., if a routine returns IINFO not equal to 0.) */
00272 
00273 /*  A       (output) COMPLEX array, dimension (LDA,max(NN)) */
00274 /*          Used to hold the matrix whose singular values are to be */
00275 /*          computed.  On exit, A contains the last matrix actually */
00276 /*          used. */
00277 
00278 /*  LDA     (input) INTEGER */
00279 /*          The leading dimension of A.  It must be at */
00280 /*          least 1 and at least max( MM ). */
00281 
00282 /*  U       (output) COMPLEX array, dimension (LDU,max(MM)) */
00283 /*          Used to hold the computed matrix of right singular vectors. */
00284 /*          On exit, U contains the last such vectors actually computed. */
00285 
00286 /*  LDU     (input) INTEGER */
00287 /*          The leading dimension of U.  It must be at */
00288 /*          least 1 and at least max( MM ). */
00289 
00290 /*  VT      (output) COMPLEX array, dimension (LDVT,max(NN)) */
00291 /*          Used to hold the computed matrix of left singular vectors. */
00292 /*          On exit, VT contains the last such vectors actually computed. */
00293 
00294 /*  LDVT    (input) INTEGER */
00295 /*          The leading dimension of VT.  It must be at */
00296 /*          least 1 and at least max( NN ). */
00297 
00298 /*  ASAV    (output) COMPLEX array, dimension (LDA,max(NN)) */
00299 /*          Used to hold a different copy of the matrix whose singular */
00300 /*          values are to be computed.  On exit, A contains the last */
00301 /*          matrix actually used. */
00302 
00303 /*  USAV    (output) COMPLEX array, dimension (LDU,max(MM)) */
00304 /*          Used to hold a different copy of the computed matrix of */
00305 /*          right singular vectors. On exit, USAV contains the last such */
00306 /*          vectors actually computed. */
00307 
00308 /*  VTSAV   (output) COMPLEX array, dimension (LDVT,max(NN)) */
00309 /*          Used to hold a different copy of the computed matrix of */
00310 /*          left singular vectors. On exit, VTSAV contains the last such */
00311 /*          vectors actually computed. */
00312 
00313 /*  S       (output) REAL array, dimension (max(min(MM,NN))) */
00314 /*          Contains the computed singular values. */
00315 
00316 /*  SSAV    (output) REAL array, dimension (max(min(MM,NN))) */
00317 /*          Contains another copy of the computed singular values. */
00318 
00319 /*  E       (output) REAL array, dimension (max(min(MM,NN))) */
00320 /*          Workspace for CGESVD. */
00321 
00322 /*  WORK    (workspace) COMPLEX array, dimension (LWORK) */
00323 
00324 /*  LWORK   (input) INTEGER */
00325 /*          The number of entries in WORK.  This must be at least */
00326 /*          MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all */
00327 /*          pairs  (M,N)=(MM(j),NN(j)) */
00328 
00329 /*  RWORK   (workspace) REAL array, */
00330 /*                      dimension ( 5*max(max(MM,NN)) ) */
00331 
00332 /*  IWORK   (workspace) INTEGER array, dimension at least 8*min(M,N) */
00333 
00334 /*  RESULT  (output) REAL array, dimension (7) */
00335 /*          The values computed by the 7 tests described above. */
00336 /*          The values are currently limited to 1/ULP, to avoid */
00337 /*          overflow. */
00338 
00339 /*  INFO    (output) INTEGER */
00340 /*          If 0, then everything ran OK. */
00341 /*           -1: NSIZES < 0 */
00342 /*           -2: Some MM(j) < 0 */
00343 /*           -3: Some NN(j) < 0 */
00344 /*           -4: NTYPES < 0 */
00345 /*           -7: THRESH < 0 */
00346 /*          -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). */
00347 /*          -12: LDU < 1 or LDU < MMAX. */
00348 /*          -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ). */
00349 /*          -21: LWORK too small. */
00350 /*          If  CLATMS, or CGESVD returns an error code, the */
00351 /*              absolute value of it is returned. */
00352 
00353 /*  ===================================================================== */
00354 
00355 /*     .. Parameters .. */
00356 /*     .. */
00357 /*     .. Local Scalars .. */
00358 /*     .. */
00359 /*     .. Local Arrays .. */
00360 /*     .. */
00361 /*     .. External Functions .. */
00362 /*     .. */
00363 /*     .. External Subroutines .. */
00364 /*     .. */
00365 /*     .. Intrinsic Functions .. */
00366 /*     .. */
00367 /*     .. Data statements .. */
00368     /* Parameter adjustments */
00369     --mm;
00370     --nn;
00371     --dotype;
00372     --iseed;
00373     asav_dim1 = *lda;
00374     asav_offset = 1 + asav_dim1;
00375     asav -= asav_offset;
00376     a_dim1 = *lda;
00377     a_offset = 1 + a_dim1;
00378     a -= a_offset;
00379     usav_dim1 = *ldu;
00380     usav_offset = 1 + usav_dim1;
00381     usav -= usav_offset;
00382     u_dim1 = *ldu;
00383     u_offset = 1 + u_dim1;
00384     u -= u_offset;
00385     vtsav_dim1 = *ldvt;
00386     vtsav_offset = 1 + vtsav_dim1;
00387     vtsav -= vtsav_offset;
00388     vt_dim1 = *ldvt;
00389     vt_offset = 1 + vt_dim1;
00390     vt -= vt_offset;
00391     --s;
00392     --ssav;
00393     --e;
00394     --work;
00395     --rwork;
00396     --iwork;
00397 
00398     /* Function Body */
00399 /*     .. */
00400 /*     .. Executable Statements .. */
00401 
00402 /*     Check for errors */
00403 
00404     *info = 0;
00405 
00406 /*     Important constants */
00407 
00408     nerrs = 0;
00409     ntestt = 0;
00410     ntestf = 0;
00411     badmm = FALSE_;
00412     badnn = FALSE_;
00413     mmax = 1;
00414     nmax = 1;
00415     mnmax = 1;
00416     minwrk = 1;
00417     i__1 = *nsizes;
00418     for (j = 1; j <= i__1; ++j) {
00419 /* Computing MAX */
00420         i__2 = mmax, i__3 = mm[j];
00421         mmax = max(i__2,i__3);
00422         if (mm[j] < 0) {
00423             badmm = TRUE_;
00424         }
00425 /* Computing MAX */
00426         i__2 = nmax, i__3 = nn[j];
00427         nmax = max(i__2,i__3);
00428         if (nn[j] < 0) {
00429             badnn = TRUE_;
00430         }
00431 /* Computing MAX */
00432 /* Computing MIN */
00433         i__4 = mm[j], i__5 = nn[j];
00434         i__2 = mnmax, i__3 = min(i__4,i__5);
00435         mnmax = max(i__2,i__3);
00436 /* Computing MAX */
00437 /* Computing MAX */
00438 /* Computing MIN */
00439         i__6 = mm[j], i__7 = nn[j];
00440 /* Computing MAX */
00441         i__9 = mm[j], i__10 = nn[j];
00442 /* Computing 2nd power */
00443         i__8 = max(i__9,i__10);
00444 /* Computing MIN */
00445         i__11 = mm[j], i__12 = nn[j];
00446 /* Computing MAX */
00447         i__13 = mm[j], i__14 = nn[j];
00448         i__4 = min(i__6,i__7) * 3 + i__8 * i__8, i__5 = min(i__11,i__12) * 5, 
00449                 i__4 = max(i__4,i__5), i__5 = max(i__13,i__14) * 3;
00450         i__2 = minwrk, i__3 = max(i__4,i__5);
00451         minwrk = max(i__2,i__3);
00452 /* L10: */
00453     }
00454 
00455 /*     Check for errors */
00456 
00457     if (*nsizes < 0) {
00458         *info = -1;
00459     } else if (badmm) {
00460         *info = -2;
00461     } else if (badnn) {
00462         *info = -3;
00463     } else if (*ntypes < 0) {
00464         *info = -4;
00465     } else if (*lda < max(1,mmax)) {
00466         *info = -10;
00467     } else if (*ldu < max(1,mmax)) {
00468         *info = -12;
00469     } else if (*ldvt < max(1,nmax)) {
00470         *info = -14;
00471     } else if (minwrk > *lwork) {
00472         *info = -21;
00473     }
00474 
00475     if (*info != 0) {
00476         i__1 = -(*info);
00477         xerbla_("CDRVBD", &i__1);
00478         return 0;
00479     }
00480 
00481 /*     Quick return if nothing to do */
00482 
00483     if (*nsizes == 0 || *ntypes == 0) {
00484         return 0;
00485     }
00486 
00487 /*     More Important constants */
00488 
00489     unfl = slamch_("S");
00490     ovfl = 1.f / unfl;
00491     ulp = slamch_("E");
00492     ulpinv = 1.f / ulp;
00493 
00494 /*     Loop over sizes, types */
00495 
00496     nerrs = 0;
00497 
00498     i__1 = *nsizes;
00499     for (jsize = 1; jsize <= i__1; ++jsize) {
00500         m = mm[jsize];
00501         n = nn[jsize];
00502         mnmin = min(m,n);
00503 
00504         if (*nsizes != 1) {
00505             mtypes = min(5,*ntypes);
00506         } else {
00507             mtypes = min(6,*ntypes);
00508         }
00509 
00510         i__2 = mtypes;
00511         for (jtype = 1; jtype <= i__2; ++jtype) {
00512             if (! dotype[jtype]) {
00513                 goto L170;
00514             }
00515             ntest = 0;
00516 
00517             for (j = 1; j <= 4; ++j) {
00518                 ioldsd[j - 1] = iseed[j];
00519 /* L20: */
00520             }
00521 
00522 /*           Compute "A" */
00523 
00524             if (mtypes > 5) {
00525                 goto L50;
00526             }
00527 
00528             if (jtype == 1) {
00529 
00530 /*              Zero matrix */
00531 
00532                 claset_("Full", &m, &n, &c_b1, &c_b1, &a[a_offset], lda);
00533                 i__3 = min(m,n);
00534                 for (i__ = 1; i__ <= i__3; ++i__) {
00535                     s[i__] = 0.f;
00536 /* L30: */
00537                 }
00538 
00539             } else if (jtype == 2) {
00540 
00541 /*              Identity matrix */
00542 
00543                 claset_("Full", &m, &n, &c_b1, &c_b2, &a[a_offset], lda);
00544                 i__3 = min(m,n);
00545                 for (i__ = 1; i__ <= i__3; ++i__) {
00546                     s[i__] = 1.f;
00547 /* L40: */
00548                 }
00549 
00550             } else {
00551 
00552 /*              (Scaled) random matrix */
00553 
00554                 if (jtype == 3) {
00555                     anorm = 1.f;
00556                 }
00557                 if (jtype == 4) {
00558                     anorm = unfl / ulp;
00559                 }
00560                 if (jtype == 5) {
00561                     anorm = ovfl * ulp;
00562                 }
00563                 r__1 = (real) mnmin;
00564                 i__3 = m - 1;
00565                 i__4 = n - 1;
00566                 clatms_(&m, &n, "U", &iseed[1], "N", &s[1], &c__4, &r__1, &
00567                         anorm, &i__3, &i__4, "N", &a[a_offset], lda, &work[1], 
00568                          &iinfo);
00569                 if (iinfo != 0) {
00570                     io___27.ciunit = *nounit;
00571                     s_wsfe(&io___27);
00572                     do_fio(&c__1, "Generator", (ftnlen)9);
00573                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00574                     do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
00575                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00576                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00577                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00578                             ;
00579                     e_wsfe();
00580                     *info = abs(iinfo);
00581                     return 0;
00582                 }
00583             }
00584 
00585 L50:
00586             clacpy_("F", &m, &n, &a[a_offset], lda, &asav[asav_offset], lda);
00587 
00588 /*           Do for minimal and adequate (for blocking) workspace */
00589 
00590             for (iwspc = 1; iwspc <= 4; ++iwspc) {
00591 
00592 /*              Test for CGESVD */
00593 
00594                 iwtmp = (min(m,n) << 1) + max(m,n);
00595                 lswork = iwtmp + (iwspc - 1) * (*lwork - iwtmp) / 3;
00596                 lswork = min(lswork,*lwork);
00597                 lswork = max(lswork,1);
00598                 if (iwspc == 4) {
00599                     lswork = *lwork;
00600                 }
00601 
00602                 for (j = 1; j <= 14; ++j) {
00603                     result[j - 1] = -1.f;
00604 /* L60: */
00605                 }
00606 
00607 /*              Factorize A */
00608 
00609                 if (iwspc > 1) {
00610                     clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
00611 , lda);
00612                 }
00613                 cgesvd_("A", "A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
00614                         usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
00615                         1], &lswork, &rwork[1], &iinfo);
00616                 if (iinfo != 0) {
00617                     io___32.ciunit = *nounit;
00618                     s_wsfe(&io___32);
00619                     do_fio(&c__1, "GESVD", (ftnlen)5);
00620                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00621                     do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
00622                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00623                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00624                     do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer));
00625                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00626                             ;
00627                     e_wsfe();
00628                     *info = abs(iinfo);
00629                     return 0;
00630                 }
00631 
00632 /*              Do tests 1--4 */
00633 
00634                 cbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
00635                         usav_offset], ldu, &ssav[1], &e[1], &vtsav[
00636                         vtsav_offset], ldvt, &work[1], &rwork[1], result);
00637                 if (m != 0 && n != 0) {
00638                     cunt01_("Columns", &mnmin, &m, &usav[usav_offset], ldu, &
00639                             work[1], lwork, &rwork[1], &result[1]);
00640                     cunt01_("Rows", &mnmin, &n, &vtsav[vtsav_offset], ldvt, &
00641                             work[1], lwork, &rwork[1], &result[2]);
00642                 }
00643                 result[3] = 0.f;
00644                 i__3 = mnmin - 1;
00645                 for (i__ = 1; i__ <= i__3; ++i__) {
00646                     if (ssav[i__] < ssav[i__ + 1]) {
00647                         result[3] = ulpinv;
00648                     }
00649                     if (ssav[i__] < 0.f) {
00650                         result[3] = ulpinv;
00651                     }
00652 /* L70: */
00653                 }
00654                 if (mnmin >= 1) {
00655                     if (ssav[mnmin] < 0.f) {
00656                         result[3] = ulpinv;
00657                     }
00658                 }
00659 
00660 /*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV */
00661 
00662                 result[4] = 0.f;
00663                 result[5] = 0.f;
00664                 result[6] = 0.f;
00665                 for (iju = 0; iju <= 3; ++iju) {
00666                     for (ijvt = 0; ijvt <= 3; ++ijvt) {
00667                         if (iju == 3 && ijvt == 3 || iju == 1 && ijvt == 1) {
00668                             goto L90;
00669                         }
00670                         *(unsigned char *)jobu = *(unsigned char *)&cjob[iju];
00671                         *(unsigned char *)jobvt = *(unsigned char *)&cjob[
00672                                 ijvt];
00673                         clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[
00674                                 a_offset], lda);
00675                         cgesvd_(jobu, jobvt, &m, &n, &a[a_offset], lda, &s[1], 
00676                                  &u[u_offset], ldu, &vt[vt_offset], ldvt, &
00677                                 work[1], &lswork, &rwork[1], &iinfo);
00678 
00679 /*                    Compare U */
00680 
00681                         dif = 0.f;
00682                         if (m > 0 && n > 0) {
00683                             if (iju == 1) {
00684                                 cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
00685                                         usav_offset], ldu, &a[a_offset], lda, 
00686                                         &work[1], lwork, &rwork[1], &dif, &
00687                                         iinfo);
00688                             } else if (iju == 2) {
00689                                 cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
00690                                         usav_offset], ldu, &u[u_offset], ldu, 
00691                                         &work[1], lwork, &rwork[1], &dif, &
00692                                         iinfo);
00693                             } else if (iju == 3) {
00694                                 cunt03_("C", &m, &m, &m, &mnmin, &usav[
00695                                         usav_offset], ldu, &u[u_offset], ldu, 
00696                                         &work[1], lwork, &rwork[1], &dif, &
00697                                         iinfo);
00698                             }
00699                         }
00700                         result[4] = dmax(result[4],dif);
00701 
00702 /*                    Compare VT */
00703 
00704                         dif = 0.f;
00705                         if (m > 0 && n > 0) {
00706                             if (ijvt == 1) {
00707                                 cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
00708                                         vtsav_offset], ldvt, &a[a_offset], 
00709                                         lda, &work[1], lwork, &rwork[1], &dif, 
00710                                          &iinfo);
00711                             } else if (ijvt == 2) {
00712                                 cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
00713                                         vtsav_offset], ldvt, &vt[vt_offset], 
00714                                         ldvt, &work[1], lwork, &rwork[1], &
00715                                         dif, &iinfo);
00716                             } else if (ijvt == 3) {
00717                                 cunt03_("R", &n, &n, &n, &mnmin, &vtsav[
00718                                         vtsav_offset], ldvt, &vt[vt_offset], 
00719                                         ldvt, &work[1], lwork, &rwork[1], &
00720                                         dif, &iinfo);
00721                             }
00722                         }
00723                         result[5] = dmax(result[5],dif);
00724 
00725 /*                    Compare S */
00726 
00727                         dif = 0.f;
00728 /* Computing MAX */
00729                         r__1 = (real) mnmin * ulp * s[1], r__2 = slamch_(
00730                                 "Safe minimum");
00731                         div = dmax(r__1,r__2);
00732                         i__3 = mnmin - 1;
00733                         for (i__ = 1; i__ <= i__3; ++i__) {
00734                             if (ssav[i__] < ssav[i__ + 1]) {
00735                                 dif = ulpinv;
00736                             }
00737                             if (ssav[i__] < 0.f) {
00738                                 dif = ulpinv;
00739                             }
00740 /* Computing MAX */
00741                             r__2 = dif, r__3 = (r__1 = ssav[i__] - s[i__], 
00742                                     dabs(r__1)) / div;
00743                             dif = dmax(r__2,r__3);
00744 /* L80: */
00745                         }
00746                         result[6] = dmax(result[6],dif);
00747 L90:
00748                         ;
00749                     }
00750 /* L100: */
00751                 }
00752 
00753 /*              Test for CGESDD */
00754 
00755                 iwtmp = (mnmin << 1) * mnmin + (mnmin << 1) + max(m,n);
00756                 lswork = iwtmp + (iwspc - 1) * (*lwork - iwtmp) / 3;
00757                 lswork = min(lswork,*lwork);
00758                 lswork = max(lswork,1);
00759                 if (iwspc == 4) {
00760                     lswork = *lwork;
00761                 }
00762 
00763 /*              Factorize A */
00764 
00765                 clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset], 
00766                         lda);
00767                 cgesdd_("A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
00768                         usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
00769                         1], &lswork, &rwork[1], &iwork[1], &iinfo);
00770                 if (iinfo != 0) {
00771                     io___39.ciunit = *nounit;
00772                     s_wsfe(&io___39);
00773                     do_fio(&c__1, "GESDD", (ftnlen)5);
00774                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00775                     do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
00776                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00777                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00778                     do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer));
00779                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00780                             ;
00781                     e_wsfe();
00782                     *info = abs(iinfo);
00783                     return 0;
00784                 }
00785 
00786 /*              Do tests 1--4 */
00787 
00788                 cbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
00789                         usav_offset], ldu, &ssav[1], &e[1], &vtsav[
00790                         vtsav_offset], ldvt, &work[1], &rwork[1], &result[7]);
00791                 if (m != 0 && n != 0) {
00792                     cunt01_("Columns", &mnmin, &m, &usav[usav_offset], ldu, &
00793                             work[1], lwork, &rwork[1], &result[8]);
00794                     cunt01_("Rows", &mnmin, &n, &vtsav[vtsav_offset], ldvt, &
00795                             work[1], lwork, &rwork[1], &result[9]);
00796                 }
00797                 result[10] = 0.f;
00798                 i__3 = mnmin - 1;
00799                 for (i__ = 1; i__ <= i__3; ++i__) {
00800                     if (ssav[i__] < ssav[i__ + 1]) {
00801                         result[10] = ulpinv;
00802                     }
00803                     if (ssav[i__] < 0.f) {
00804                         result[10] = ulpinv;
00805                     }
00806 /* L110: */
00807                 }
00808                 if (mnmin >= 1) {
00809                     if (ssav[mnmin] < 0.f) {
00810                         result[10] = ulpinv;
00811                     }
00812                 }
00813 
00814 /*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV */
00815 
00816                 result[11] = 0.f;
00817                 result[12] = 0.f;
00818                 result[13] = 0.f;
00819                 for (ijq = 0; ijq <= 2; ++ijq) {
00820                     *(unsigned char *)jobq = *(unsigned char *)&cjob[ijq];
00821                     clacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
00822 , lda);
00823                     cgesdd_(jobq, &m, &n, &a[a_offset], lda, &s[1], &u[
00824                             u_offset], ldu, &vt[vt_offset], ldvt, &work[1], &
00825                             lswork, &rwork[1], &iwork[1], &iinfo);
00826 
00827 /*                 Compare U */
00828 
00829                     dif = 0.f;
00830                     if (m > 0 && n > 0) {
00831                         if (ijq == 1) {
00832                             if (m >= n) {
00833                                 cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
00834                                         usav_offset], ldu, &a[a_offset], lda, 
00835                                         &work[1], lwork, &rwork[1], &dif, &
00836                                         iinfo);
00837                             } else {
00838                                 cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
00839                                         usav_offset], ldu, &u[u_offset], ldu, 
00840                                         &work[1], lwork, &rwork[1], &dif, &
00841                                         iinfo);
00842                             }
00843                         } else if (ijq == 2) {
00844                             cunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
00845                                     usav_offset], ldu, &u[u_offset], ldu, &
00846                                     work[1], lwork, &rwork[1], &dif, &iinfo);
00847                         }
00848                     }
00849                     result[11] = dmax(result[11],dif);
00850 
00851 /*                 Compare VT */
00852 
00853                     dif = 0.f;
00854                     if (m > 0 && n > 0) {
00855                         if (ijq == 1) {
00856                             if (m >= n) {
00857                                 cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
00858                                         vtsav_offset], ldvt, &vt[vt_offset], 
00859                                         ldvt, &work[1], lwork, &rwork[1], &
00860                                         dif, &iinfo);
00861                             } else {
00862                                 cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
00863                                         vtsav_offset], ldvt, &a[a_offset], 
00864                                         lda, &work[1], lwork, &rwork[1], &dif, 
00865                                          &iinfo);
00866                             }
00867                         } else if (ijq == 2) {
00868                             cunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
00869                                     vtsav_offset], ldvt, &vt[vt_offset], ldvt, 
00870                                      &work[1], lwork, &rwork[1], &dif, &iinfo);
00871                         }
00872                     }
00873                     result[12] = dmax(result[12],dif);
00874 
00875 /*                 Compare S */
00876 
00877                     dif = 0.f;
00878 /* Computing MAX */
00879                     r__1 = (real) mnmin * ulp * s[1], r__2 = slamch_("Safe m"
00880                             "inimum");
00881                     div = dmax(r__1,r__2);
00882                     i__3 = mnmin - 1;
00883                     for (i__ = 1; i__ <= i__3; ++i__) {
00884                         if (ssav[i__] < ssav[i__ + 1]) {
00885                             dif = ulpinv;
00886                         }
00887                         if (ssav[i__] < 0.f) {
00888                             dif = ulpinv;
00889                         }
00890 /* Computing MAX */
00891                         r__2 = dif, r__3 = (r__1 = ssav[i__] - s[i__], dabs(
00892                                 r__1)) / div;
00893                         dif = dmax(r__2,r__3);
00894 /* L120: */
00895                     }
00896                     result[13] = dmax(result[13],dif);
00897 /* L130: */
00898                 }
00899 
00900 /*              End of Loop -- Check for RESULT(j) > THRESH */
00901 
00902                 ntest = 0;
00903                 nfail = 0;
00904                 for (j = 1; j <= 14; ++j) {
00905                     if (result[j - 1] >= 0.f) {
00906                         ++ntest;
00907                     }
00908                     if (result[j - 1] >= *thresh) {
00909                         ++nfail;
00910                     }
00911 /* L140: */
00912                 }
00913 
00914                 if (nfail > 0) {
00915                     ++ntestf;
00916                 }
00917                 if (ntestf == 1) {
00918                     io___43.ciunit = *nounit;
00919                     s_wsfe(&io___43);
00920                     e_wsfe();
00921                     io___44.ciunit = *nounit;
00922                     s_wsfe(&io___44);
00923                     do_fio(&c__1, (char *)&(*thresh), (ftnlen)sizeof(real));
00924                     e_wsfe();
00925                     ntestf = 2;
00926                 }
00927 
00928                 for (j = 1; j <= 14; ++j) {
00929                     if (result[j - 1] >= *thresh) {
00930                         io___45.ciunit = *nounit;
00931                         s_wsfe(&io___45);
00932                         do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
00933                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00934                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
00935                                 ;
00936                         do_fio(&c__1, (char *)&iwspc, (ftnlen)sizeof(integer))
00937                                 ;
00938                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
00939                                 integer));
00940                         do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
00941                         do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(
00942                                 real));
00943                         e_wsfe();
00944                     }
00945 /* L150: */
00946                 }
00947 
00948                 nerrs += nfail;
00949                 ntestt += ntest;
00950 
00951 /* L160: */
00952             }
00953 
00954 L170:
00955             ;
00956         }
00957 /* L180: */
00958     }
00959 
00960 /*     Summary */
00961 
00962     alasvm_("CBD", nounit, &nerrs, &ntestt, &c__0);
00963 
00964 
00965     return 0;
00966 
00967 /*     End of CDRVBD */
00968 
00969 } /* cdrvbd_ */


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