zchkbb.c
Go to the documentation of this file.
00001 /* zchkbb.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 doublecomplex c_b1 = {0.,0.};
00019 static doublecomplex c_b2 = {1.,0.};
00020 static integer c__0 = 0;
00021 static integer c__6 = 6;
00022 static doublereal c_b33 = 1.;
00023 static integer c__1 = 1;
00024 static doublereal c_b41 = 0.;
00025 static integer c__4 = 4;
00026 static integer c_n1 = -1;
00027 
00028 /* Subroutine */ int zchkbb_(integer *nsizes, integer *mval, integer *nval, 
00029         integer *nwdths, integer *kk, integer *ntypes, logical *dotype, 
00030         integer *nrhs, integer *iseed, doublereal *thresh, integer *nounit, 
00031         doublecomplex *a, integer *lda, doublecomplex *ab, integer *ldab, 
00032         doublereal *bd, doublereal *be, doublecomplex *q, integer *ldq, 
00033         doublecomplex *p, integer *ldp, doublecomplex *c__, integer *ldc, 
00034         doublecomplex *cc, doublecomplex *work, integer *lwork, doublereal *
00035         rwork, doublereal *result, integer *info)
00036 {
00037     /* Initialized data */
00038 
00039     static integer ktype[15] = { 1,2,4,4,4,4,4,6,6,6,6,6,9,9,9 };
00040     static integer kmagn[15] = { 1,1,1,1,1,2,3,1,1,1,2,3,1,2,3 };
00041     static integer kmode[15] = { 0,0,4,3,1,4,4,4,3,1,4,4,0,0,0 };
00042 
00043     /* Format strings */
00044     static char fmt_9999[] = "(\002 ZCHKBB: \002,a,\002 returned INFO=\002,i"
00045             "5,\002.\002,/9x,\002M=\002,i5,\002 N=\002,i5,\002 K=\002,i5,\002"
00046             ", JTYPE=\002,i5,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
00047     static char fmt_9998[] = "(\002 M =\002,i4,\002 N=\002,i4,\002, K=\002,i"
00048             "3,\002, seed=\002,4(i4,\002,\002),\002 type \002,i2,\002, test"
00049             "(\002,i2,\002)=\002,g10.3)";
00050 
00051     /* System generated locals */
00052     integer a_dim1, a_offset, ab_dim1, ab_offset, c_dim1, c_offset, cc_dim1, 
00053             cc_offset, p_dim1, p_offset, q_dim1, q_offset, i__1, i__2, i__3, 
00054             i__4, i__5, i__6, i__7, i__8, i__9;
00055 
00056     /* Builtin functions */
00057     double sqrt(doublereal);
00058     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00059 
00060     /* Local variables */
00061     integer i__, j, k, m, n, kl, jr, ku;
00062     doublereal ulp, cond;
00063     integer jcol, kmax, mmax, nmax;
00064     doublereal unfl, ovfl;
00065     logical badmm, badnn;
00066     integer imode, iinfo;
00067     extern /* Subroutine */ int zbdt01_(integer *, integer *, integer *, 
00068             doublecomplex *, integer *, doublecomplex *, integer *, 
00069             doublereal *, doublereal *, doublecomplex *, integer *, 
00070             doublecomplex *, doublereal *, doublereal *), zbdt02_(integer *, 
00071             integer *, doublecomplex *, integer *, doublecomplex *, integer *, 
00072              doublecomplex *, integer *, doublecomplex *, doublereal *, 
00073             doublereal *);
00074     doublereal anorm;
00075     integer mnmin, mnmax, nmats, jsize, nerrs, itype, jtype, ntest;
00076     extern /* Subroutine */ int dlahd2_(integer *, char *), zunt01_(
00077             char *, integer *, integer *, doublecomplex *, integer *, 
00078             doublecomplex *, integer *, doublereal *, doublereal *);
00079     logical badnnb;
00080     extern doublereal dlamch_(char *);
00081     extern /* Subroutine */ int zgbbrd_(char *, integer *, integer *, integer 
00082             *, integer *, integer *, doublecomplex *, integer *, doublereal *, 
00083              doublereal *, doublecomplex *, integer *, doublecomplex *, 
00084             integer *, doublecomplex *, integer *, doublecomplex *, 
00085             doublereal *, integer *);
00086     integer idumma[1];
00087     extern /* Subroutine */ int xerbla_(char *, integer *);
00088     integer ioldsd[4];
00089     extern /* Subroutine */ int dlasum_(char *, integer *, integer *, integer 
00090             *);
00091     doublereal amninv;
00092     integer jwidth;
00093     extern /* Subroutine */ int zlacpy_(char *, integer *, integer *, 
00094             doublecomplex *, integer *, doublecomplex *, integer *), 
00095             zlaset_(char *, integer *, integer *, doublecomplex *, 
00096             doublecomplex *, doublecomplex *, integer *), zlatmr_(
00097             integer *, integer *, char *, integer *, char *, doublecomplex *, 
00098             integer *, doublereal *, doublecomplex *, char *, char *, 
00099             doublecomplex *, integer *, doublereal *, doublecomplex *, 
00100             integer *, doublereal *, char *, integer *, integer *, integer *, 
00101             doublereal *, doublereal *, char *, doublecomplex *, integer *, 
00102             integer *, integer *);
00103     doublereal rtunfl, rtovfl, ulpinv;
00104     extern /* Subroutine */ int zlatms_(integer *, integer *, char *, integer 
00105             *, char *, doublereal *, integer *, doublereal *, doublereal *, 
00106             integer *, integer *, char *, doublecomplex *, integer *, 
00107             doublecomplex *, integer *);
00108     integer mtypes, ntestt;
00109 
00110     /* Fortran I/O blocks */
00111     static cilist io___41 = { 0, 0, 0, fmt_9999, 0 };
00112     static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };
00113     static cilist io___45 = { 0, 0, 0, fmt_9998, 0 };
00114 
00115 
00116 
00117 /*  -- LAPACK test routine (new routine for release 2.0) -- */
00118 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00119 /*     November 2006 */
00120 
00121 /*     .. Scalar Arguments .. */
00122 /*     .. */
00123 /*     .. Array Arguments .. */
00124 /*     .. */
00125 
00126 /*  Purpose */
00127 /*  ======= */
00128 
00129 /*  ZCHKBB tests the reduction of a general complex rectangular band */
00130 /*  matrix to real bidiagonal form. */
00131 
00132 /*  ZGBBRD factors a general band matrix A as  Q B P* , where * means */
00133 /*  conjugate transpose, B is upper bidiagonal, and Q and P are unitary; */
00134 /*  ZGBBRD can also overwrite a given matrix C with Q* C . */
00135 
00136 /*  For each pair of matrix dimensions (M,N) and each selected matrix */
00137 /*  type, an M by N matrix A and an M by NRHS matrix C are generated. */
00138 /*  The problem dimensions are as follows */
00139 /*     A:          M x N */
00140 /*     Q:          M x M */
00141 /*     P:          N x N */
00142 /*     B:          min(M,N) x min(M,N) */
00143 /*     C:          M x NRHS */
00144 
00145 /*  For each generated matrix, 4 tests are performed: */
00146 
00147 /*  (1)   | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P' */
00148 
00149 /*  (2)   | I - Q' Q | / ( M ulp ) */
00150 
00151 /*  (3)   | I - PT PT' | / ( N ulp ) */
00152 
00153 /*  (4)   | Y - Q' C | / ( |Y| max(M,NRHS) ulp ), where Y = Q' C. */
00154 
00155 /*  The "types" are specified by a logical array DOTYPE( 1:NTYPES ); */
00156 /*  if DOTYPE(j) is .TRUE., then matrix type "j" will be generated. */
00157 /*  Currently, the list of possible types is: */
00158 
00159 /*  The possible matrix types are */
00160 
00161 /*  (1)  The zero matrix. */
00162 /*  (2)  The identity matrix. */
00163 
00164 /*  (3)  A diagonal matrix with evenly spaced entries */
00165 /*       1, ..., ULP  and random signs. */
00166 /*       (ULP = (first number larger than 1) - 1 ) */
00167 /*  (4)  A diagonal matrix with geometrically spaced entries */
00168 /*       1, ..., ULP  and random signs. */
00169 /*  (5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP */
00170 /*       and random signs. */
00171 
00172 /*  (6)  Same as (3), but multiplied by SQRT( overflow threshold ) */
00173 /*  (7)  Same as (3), but multiplied by SQRT( underflow threshold ) */
00174 
00175 /*  (8)  A matrix of the form  U D V, where U and V are orthogonal and */
00176 /*       D has evenly spaced entries 1, ..., ULP with random signs */
00177 /*       on the diagonal. */
00178 
00179 /*  (9)  A matrix of the form  U D V, where U and V are orthogonal and */
00180 /*       D has geometrically spaced entries 1, ..., ULP with random */
00181 /*       signs on the diagonal. */
00182 
00183 /*  (10) A matrix of the form  U D V, where U and V are orthogonal and */
00184 /*       D has "clustered" entries 1, ULP,..., ULP with random */
00185 /*       signs on the diagonal. */
00186 
00187 /*  (11) Same as (8), but multiplied by SQRT( overflow threshold ) */
00188 /*  (12) Same as (8), but multiplied by SQRT( underflow threshold ) */
00189 
00190 /*  (13) Rectangular matrix with random entries chosen from (-1,1). */
00191 /*  (14) Same as (13), but multiplied by SQRT( overflow threshold ) */
00192 /*  (15) Same as (13), but multiplied by SQRT( underflow threshold ) */
00193 
00194 /*  Arguments */
00195 /*  ========= */
00196 
00197 /*  NSIZES  (input) INTEGER */
00198 /*          The number of values of M and N contained in the vectors */
00199 /*          MVAL and NVAL.  The matrix sizes are used in pairs (M,N). */
00200 /*          If NSIZES is zero, ZCHKBB does nothing.  NSIZES must be at */
00201 /*          least zero. */
00202 
00203 /*  MVAL    (input) INTEGER array, dimension (NSIZES) */
00204 /*          The values of the matrix row dimension M. */
00205 
00206 /*  NVAL    (input) INTEGER array, dimension (NSIZES) */
00207 /*          The values of the matrix column dimension N. */
00208 
00209 /*  NWDTHS  (input) INTEGER */
00210 /*          The number of bandwidths to use.  If it is zero, */
00211 /*          ZCHKBB does nothing.  It must be at least zero. */
00212 
00213 /*  KK      (input) INTEGER array, dimension (NWDTHS) */
00214 /*          An array containing the bandwidths to be used for the band */
00215 /*          matrices.  The values must be at least zero. */
00216 
00217 /*  NTYPES  (input) INTEGER */
00218 /*          The number of elements in DOTYPE.   If it is zero, ZCHKBB */
00219 /*          does nothing.  It must be at least zero.  If it is MAXTYP+1 */
00220 /*          and NSIZES is 1, then an additional type, MAXTYP+1 is */
00221 /*          defined, which is to use whatever matrix is in A.  This */
00222 /*          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and */
00223 /*          DOTYPE(MAXTYP+1) is .TRUE. . */
00224 
00225 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00226 /*          If DOTYPE(j) is .TRUE., then for each size in NN a */
00227 /*          matrix of that size and of type j will be generated. */
00228 /*          If NTYPES is smaller than the maximum number of types */
00229 /*          defined (PARAMETER MAXTYP), then types NTYPES+1 through */
00230 /*          MAXTYP will not be generated.  If NTYPES is larger */
00231 /*          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) */
00232 /*          will be ignored. */
00233 
00234 /*  NRHS    (input) INTEGER */
00235 /*          The number of columns in the "right-hand side" matrix C. */
00236 /*          If NRHS = 0, then the operations on the right-hand side will */
00237 /*          not be tested. NRHS must be at least 0. */
00238 
00239 /*  ISEED   (input/output) INTEGER array, dimension (4) */
00240 /*          On entry ISEED specifies the seed of the random number */
00241 /*          generator. The array elements should be between 0 and 4095; */
00242 /*          if not they will be reduced mod 4096.  Also, ISEED(4) must */
00243 /*          be odd.  The random number generator uses a linear */
00244 /*          congruential sequence limited to small integers, and so */
00245 /*          should produce machine independent random numbers. The */
00246 /*          values of ISEED are changed on exit, and can be used in the */
00247 /*          next call to ZCHKBB to continue the same random number */
00248 /*          sequence. */
00249 
00250 /*  THRESH  (input) DOUBLE PRECISION */
00251 /*          A test will count as "failed" if the "error", computed as */
00252 /*          described above, exceeds THRESH.  Note that the error */
00253 /*          is scaled to be O(1), so THRESH should be a reasonably */
00254 /*          small multiple of 1, e.g., 10 or 100.  In particular, */
00255 /*          it should not depend on the precision (single vs. double) */
00256 /*          or the size of the matrix.  It must be at least zero. */
00257 
00258 /*  NOUNIT  (input) INTEGER */
00259 /*          The FORTRAN unit number for printing out error messages */
00260 /*          (e.g., if a routine returns IINFO not equal to 0.) */
00261 
00262 /*  A       (input/workspace) DOUBLE PRECISION array, dimension */
00263 /*                            (LDA, max(NN)) */
00264 /*          Used to hold the matrix A. */
00265 
00266 /*  LDA     (input) INTEGER */
00267 /*          The leading dimension of A.  It must be at least 1 */
00268 /*          and at least max( NN ). */
00269 
00270 /*  AB      (workspace) DOUBLE PRECISION array, dimension (LDAB, max(NN)) */
00271 /*          Used to hold A in band storage format. */
00272 
00273 /*  LDAB    (input) INTEGER */
00274 /*          The leading dimension of AB.  It must be at least 2 (not 1!) */
00275 /*          and at least max( KK )+1. */
00276 
00277 /*  BD      (workspace) DOUBLE PRECISION array, dimension (max(NN)) */
00278 /*          Used to hold the diagonal of the bidiagonal matrix computed */
00279 /*          by ZGBBRD. */
00280 
00281 /*  BE      (workspace) DOUBLE PRECISION array, dimension (max(NN)) */
00282 /*          Used to hold the off-diagonal of the bidiagonal matrix */
00283 /*          computed by ZGBBRD. */
00284 
00285 /*  Q       (workspace) COMPLEX*16 array, dimension (LDQ, max(NN)) */
00286 /*          Used to hold the unitary matrix Q computed by ZGBBRD. */
00287 
00288 /*  LDQ     (input) INTEGER */
00289 /*          The leading dimension of Q.  It must be at least 1 */
00290 /*          and at least max( NN ). */
00291 
00292 /*  P       (workspace) COMPLEX*16 array, dimension (LDP, max(NN)) */
00293 /*          Used to hold the unitary matrix P computed by ZGBBRD. */
00294 
00295 /*  LDP     (input) INTEGER */
00296 /*          The leading dimension of P.  It must be at least 1 */
00297 /*          and at least max( NN ). */
00298 
00299 /*  C       (workspace) COMPLEX*16 array, dimension (LDC, max(NN)) */
00300 /*          Used to hold the matrix C updated by ZGBBRD. */
00301 
00302 /*  LDC     (input) INTEGER */
00303 /*          The leading dimension of U.  It must be at least 1 */
00304 /*          and at least max( NN ). */
00305 
00306 /*  CC      (workspace) COMPLEX*16 array, dimension (LDC, max(NN)) */
00307 /*          Used to hold a copy of the matrix C. */
00308 
00309 /*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK) */
00310 
00311 /*  LWORK   (input) INTEGER */
00312 /*          The number of entries in WORK.  This must be at least */
00313 /*          max( LDA+1, max(NN)+1 )*max(NN). */
00314 
00315 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (max(NN)) */
00316 
00317 /*  RESULT  (output) DOUBLE PRECISION array, dimension (4) */
00318 /*          The values computed by the tests described above. */
00319 /*          The values are currently limited to 1/ulp, to avoid */
00320 /*          overflow. */
00321 
00322 /*  INFO    (output) INTEGER */
00323 /*          If 0, then everything ran OK. */
00324 
00325 /* ----------------------------------------------------------------------- */
00326 
00327 /*       Some Local Variables and Parameters: */
00328 /*       ---- ----- --------- --- ---------- */
00329 /*       ZERO, ONE       Real 0 and 1. */
00330 /*       MAXTYP          The number of types defined. */
00331 /*       NTEST           The number of tests performed, or which can */
00332 /*                       be performed so far, for the current matrix. */
00333 /*       NTESTT          The total number of tests performed so far. */
00334 /*       NMAX            Largest value in NN. */
00335 /*       NMATS           The number of matrices generated so far. */
00336 /*       NERRS           The number of tests which have exceeded THRESH */
00337 /*                       so far. */
00338 /*       COND, IMODE     Values to be passed to the matrix generators. */
00339 /*       ANORM           Norm of A; passed to matrix generators. */
00340 
00341 /*       OVFL, UNFL      Overflow and underflow thresholds. */
00342 /*       ULP, ULPINV     Finest relative precision and its inverse. */
00343 /*       RTOVFL, RTUNFL  Square roots of the previous 2 values. */
00344 /*               The following four arrays decode JTYPE: */
00345 /*       KTYPE(j)        The general type (1-10) for type "j". */
00346 /*       KMODE(j)        The MODE value to be passed to the matrix */
00347 /*                       generator for type "j". */
00348 /*       KMAGN(j)        The order of magnitude ( O(1), */
00349 /*                       O(overflow^(1/2) ), O(underflow^(1/2) ) */
00350 
00351 /*  ===================================================================== */
00352 
00353 /*     .. Parameters .. */
00354 /*     .. */
00355 /*     .. Local Scalars .. */
00356 /*     .. */
00357 /*     .. Local Arrays .. */
00358 /*     .. */
00359 /*     .. External Functions .. */
00360 /*     .. */
00361 /*     .. External Subroutines .. */
00362 /*     .. */
00363 /*     .. Intrinsic Functions .. */
00364 /*     .. */
00365 /*     .. Data statements .. */
00366     /* Parameter adjustments */
00367     --mval;
00368     --nval;
00369     --kk;
00370     --dotype;
00371     --iseed;
00372     a_dim1 = *lda;
00373     a_offset = 1 + a_dim1;
00374     a -= a_offset;
00375     ab_dim1 = *ldab;
00376     ab_offset = 1 + ab_dim1;
00377     ab -= ab_offset;
00378     --bd;
00379     --be;
00380     q_dim1 = *ldq;
00381     q_offset = 1 + q_dim1;
00382     q -= q_offset;
00383     p_dim1 = *ldp;
00384     p_offset = 1 + p_dim1;
00385     p -= p_offset;
00386     cc_dim1 = *ldc;
00387     cc_offset = 1 + cc_dim1;
00388     cc -= cc_offset;
00389     c_dim1 = *ldc;
00390     c_offset = 1 + c_dim1;
00391     c__ -= c_offset;
00392     --work;
00393     --rwork;
00394     --result;
00395 
00396     /* Function Body */
00397 /*     .. */
00398 /*     .. Executable Statements .. */
00399 
00400 /*     Check for errors */
00401 
00402     ntestt = 0;
00403     *info = 0;
00404 
00405 /*     Important constants */
00406 
00407     badmm = FALSE_;
00408     badnn = FALSE_;
00409     mmax = 1;
00410     nmax = 1;
00411     mnmax = 1;
00412     i__1 = *nsizes;
00413     for (j = 1; j <= i__1; ++j) {
00414 /* Computing MAX */
00415         i__2 = mmax, i__3 = mval[j];
00416         mmax = max(i__2,i__3);
00417         if (mval[j] < 0) {
00418             badmm = TRUE_;
00419         }
00420 /* Computing MAX */
00421         i__2 = nmax, i__3 = nval[j];
00422         nmax = max(i__2,i__3);
00423         if (nval[j] < 0) {
00424             badnn = TRUE_;
00425         }
00426 /* Computing MAX */
00427 /* Computing MIN */
00428         i__4 = mval[j], i__5 = nval[j];
00429         i__2 = mnmax, i__3 = min(i__4,i__5);
00430         mnmax = max(i__2,i__3);
00431 /* L10: */
00432     }
00433 
00434     badnnb = FALSE_;
00435     kmax = 0;
00436     i__1 = *nwdths;
00437     for (j = 1; j <= i__1; ++j) {
00438 /* Computing MAX */
00439         i__2 = kmax, i__3 = kk[j];
00440         kmax = max(i__2,i__3);
00441         if (kk[j] < 0) {
00442             badnnb = TRUE_;
00443         }
00444 /* L20: */
00445     }
00446 
00447 /*     Check for errors */
00448 
00449     if (*nsizes < 0) {
00450         *info = -1;
00451     } else if (badmm) {
00452         *info = -2;
00453     } else if (badnn) {
00454         *info = -3;
00455     } else if (*nwdths < 0) {
00456         *info = -4;
00457     } else if (badnnb) {
00458         *info = -5;
00459     } else if (*ntypes < 0) {
00460         *info = -6;
00461     } else if (*nrhs < 0) {
00462         *info = -8;
00463     } else if (*lda < nmax) {
00464         *info = -13;
00465     } else if (*ldab < (kmax << 1) + 1) {
00466         *info = -15;
00467     } else if (*ldq < nmax) {
00468         *info = -19;
00469     } else if (*ldp < nmax) {
00470         *info = -21;
00471     } else if (*ldc < nmax) {
00472         *info = -23;
00473     } else if ((max(*lda,nmax) + 1) * nmax > *lwork) {
00474         *info = -26;
00475     }
00476 
00477     if (*info != 0) {
00478         i__1 = -(*info);
00479         xerbla_("ZCHKBB", &i__1);
00480         return 0;
00481     }
00482 
00483 /*     Quick return if possible */
00484 
00485     if (*nsizes == 0 || *ntypes == 0 || *nwdths == 0) {
00486         return 0;
00487     }
00488 
00489 /*     More Important constants */
00490 
00491     unfl = dlamch_("Safe minimum");
00492     ovfl = 1. / unfl;
00493     ulp = dlamch_("Epsilon") * dlamch_("Base");
00494     ulpinv = 1. / ulp;
00495     rtunfl = sqrt(unfl);
00496     rtovfl = sqrt(ovfl);
00497 
00498 /*     Loop over sizes, widths, types */
00499 
00500     nerrs = 0;
00501     nmats = 0;
00502 
00503     i__1 = *nsizes;
00504     for (jsize = 1; jsize <= i__1; ++jsize) {
00505         m = mval[jsize];
00506         n = nval[jsize];
00507         mnmin = min(m,n);
00508 /* Computing MAX */
00509         i__2 = max(1,m);
00510         amninv = 1. / (doublereal) max(i__2,n);
00511 
00512         i__2 = *nwdths;
00513         for (jwidth = 1; jwidth <= i__2; ++jwidth) {
00514             k = kk[jwidth];
00515             if (k >= m && k >= n) {
00516                 goto L150;
00517             }
00518 /* Computing MAX */
00519 /* Computing MIN */
00520             i__5 = m - 1;
00521             i__3 = 0, i__4 = min(i__5,k);
00522             kl = max(i__3,i__4);
00523 /* Computing MAX */
00524 /* Computing MIN */
00525             i__5 = n - 1;
00526             i__3 = 0, i__4 = min(i__5,k);
00527             ku = max(i__3,i__4);
00528 
00529             if (*nsizes != 1) {
00530                 mtypes = min(15,*ntypes);
00531             } else {
00532                 mtypes = min(16,*ntypes);
00533             }
00534 
00535             i__3 = mtypes;
00536             for (jtype = 1; jtype <= i__3; ++jtype) {
00537                 if (! dotype[jtype]) {
00538                     goto L140;
00539                 }
00540                 ++nmats;
00541                 ntest = 0;
00542 
00543                 for (j = 1; j <= 4; ++j) {
00544                     ioldsd[j - 1] = iseed[j];
00545 /* L30: */
00546                 }
00547 
00548 /*              Compute "A". */
00549 
00550 /*              Control parameters: */
00551 
00552 /*                  KMAGN  KMODE        KTYPE */
00553 /*              =1  O(1)   clustered 1  zero */
00554 /*              =2  large  clustered 2  identity */
00555 /*              =3  small  exponential  (none) */
00556 /*              =4         arithmetic   diagonal, (w/ singular values) */
00557 /*              =5         random log   (none) */
00558 /*              =6         random       nonhermitian, w/ singular values */
00559 /*              =7                      (none) */
00560 /*              =8                      (none) */
00561 /*              =9                      random nonhermitian */
00562 
00563                 if (mtypes > 15) {
00564                     goto L90;
00565                 }
00566 
00567                 itype = ktype[jtype - 1];
00568                 imode = kmode[jtype - 1];
00569 
00570 /*              Compute norm */
00571 
00572                 switch (kmagn[jtype - 1]) {
00573                     case 1:  goto L40;
00574                     case 2:  goto L50;
00575                     case 3:  goto L60;
00576                 }
00577 
00578 L40:
00579                 anorm = 1.;
00580                 goto L70;
00581 
00582 L50:
00583                 anorm = rtovfl * ulp * amninv;
00584                 goto L70;
00585 
00586 L60:
00587                 anorm = rtunfl * max(m,n) * ulpinv;
00588                 goto L70;
00589 
00590 L70:
00591 
00592                 zlaset_("Full", lda, &n, &c_b1, &c_b1, &a[a_offset], lda);
00593                 zlaset_("Full", ldab, &n, &c_b1, &c_b1, &ab[ab_offset], ldab);
00594                 iinfo = 0;
00595                 cond = ulpinv;
00596 
00597 /*              Special Matrices -- Identity & Jordan block */
00598 
00599 /*                 Zero */
00600 
00601                 if (itype == 1) {
00602                     iinfo = 0;
00603 
00604                 } else if (itype == 2) {
00605 
00606 /*                 Identity */
00607 
00608                     i__4 = n;
00609                     for (jcol = 1; jcol <= i__4; ++jcol) {
00610                         i__5 = jcol + jcol * a_dim1;
00611                         a[i__5].r = anorm, a[i__5].i = 0.;
00612 /* L80: */
00613                     }
00614 
00615                 } else if (itype == 4) {
00616 
00617 /*                 Diagonal Matrix, singular values specified */
00618 
00619                     zlatms_(&m, &n, "S", &iseed[1], "N", &rwork[1], &imode, &
00620                             cond, &anorm, &c__0, &c__0, "N", &a[a_offset], 
00621                             lda, &work[1], &iinfo);
00622 
00623                 } else if (itype == 6) {
00624 
00625 /*                 Nonhermitian, singular values specified */
00626 
00627                     zlatms_(&m, &n, "S", &iseed[1], "N", &rwork[1], &imode, &
00628                             cond, &anorm, &kl, &ku, "N", &a[a_offset], lda, &
00629                             work[1], &iinfo);
00630 
00631                 } else if (itype == 9) {
00632 
00633 /*                 Nonhermitian, random entries */
00634 
00635                     zlatmr_(&m, &n, "S", &iseed[1], "N", &work[1], &c__6, &
00636                             c_b33, &c_b2, "T", "N", &work[n + 1], &c__1, &
00637                             c_b33, &work[(n << 1) + 1], &c__1, &c_b33, "N", 
00638                             idumma, &kl, &ku, &c_b41, &anorm, "N", &a[
00639                             a_offset], lda, idumma, &iinfo);
00640 
00641                 } else {
00642 
00643                     iinfo = 1;
00644                 }
00645 
00646 /*              Generate Right-Hand Side */
00647 
00648                 zlatmr_(&m, nrhs, "S", &iseed[1], "N", &work[1], &c__6, &
00649                         c_b33, &c_b2, "T", "N", &work[m + 1], &c__1, &c_b33, &
00650                         work[(m << 1) + 1], &c__1, &c_b33, "N", idumma, &m, 
00651                         nrhs, &c_b41, &c_b33, "NO", &c__[c_offset], ldc, 
00652                         idumma, &iinfo);
00653 
00654                 if (iinfo != 0) {
00655                     io___41.ciunit = *nounit;
00656                     s_wsfe(&io___41);
00657                     do_fio(&c__1, "Generator", (ftnlen)9);
00658                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00659                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00660                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00661                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00662                             ;
00663                     e_wsfe();
00664                     *info = abs(iinfo);
00665                     return 0;
00666                 }
00667 
00668 L90:
00669 
00670 /*              Copy A to band storage. */
00671 
00672                 i__4 = n;
00673                 for (j = 1; j <= i__4; ++j) {
00674 /* Computing MAX */
00675                     i__5 = 1, i__6 = j - ku;
00676 /* Computing MIN */
00677                     i__8 = m, i__9 = j + kl;
00678                     i__7 = min(i__8,i__9);
00679                     for (i__ = max(i__5,i__6); i__ <= i__7; ++i__) {
00680                         i__5 = ku + 1 + i__ - j + j * ab_dim1;
00681                         i__6 = i__ + j * a_dim1;
00682                         ab[i__5].r = a[i__6].r, ab[i__5].i = a[i__6].i;
00683 /* L100: */
00684                     }
00685 /* L110: */
00686                 }
00687 
00688 /*              Copy C */
00689 
00690                 zlacpy_("Full", &m, nrhs, &c__[c_offset], ldc, &cc[cc_offset], 
00691                          ldc);
00692 
00693 /*              Call ZGBBRD to compute B, Q and P, and to update C. */
00694 
00695                 zgbbrd_("B", &m, &n, nrhs, &kl, &ku, &ab[ab_offset], ldab, &
00696                         bd[1], &be[1], &q[q_offset], ldq, &p[p_offset], ldp, &
00697                         cc[cc_offset], ldc, &work[1], &rwork[1], &iinfo);
00698 
00699                 if (iinfo != 0) {
00700                     io___43.ciunit = *nounit;
00701                     s_wsfe(&io___43);
00702                     do_fio(&c__1, "ZGBBRD", (ftnlen)6);
00703                     do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
00704                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00705                     do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
00706                     do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
00707                             ;
00708                     e_wsfe();
00709                     *info = abs(iinfo);
00710                     if (iinfo < 0) {
00711                         return 0;
00712                     } else {
00713                         result[1] = ulpinv;
00714                         goto L120;
00715                     }
00716                 }
00717 
00718 /*              Test 1:  Check the decomposition A := Q * B * P' */
00719 /*                   2:  Check the orthogonality of Q */
00720 /*                   3:  Check the orthogonality of P */
00721 /*                   4:  Check the computation of Q' * C */
00722 
00723                 zbdt01_(&m, &n, &c_n1, &a[a_offset], lda, &q[q_offset], ldq, &
00724                         bd[1], &be[1], &p[p_offset], ldp, &work[1], &rwork[1], 
00725                          &result[1]);
00726                 zunt01_("Columns", &m, &m, &q[q_offset], ldq, &work[1], lwork, 
00727                          &rwork[1], &result[2]);
00728                 zunt01_("Rows", &n, &n, &p[p_offset], ldp, &work[1], lwork, &
00729                         rwork[1], &result[3]);
00730                 zbdt02_(&m, nrhs, &c__[c_offset], ldc, &cc[cc_offset], ldc, &
00731                         q[q_offset], ldq, &work[1], &rwork[1], &result[4]);
00732 
00733 /*              End of Loop -- Check for RESULT(j) > THRESH */
00734 
00735                 ntest = 4;
00736 L120:
00737                 ntestt += ntest;
00738 
00739 /*              Print out tests which fail. */
00740 
00741                 i__4 = ntest;
00742                 for (jr = 1; jr <= i__4; ++jr) {
00743                     if (result[jr] >= *thresh) {
00744                         if (nerrs == 0) {
00745                             dlahd2_(nounit, "ZBB");
00746                         }
00747                         ++nerrs;
00748                         io___45.ciunit = *nounit;
00749                         s_wsfe(&io___45);
00750                         do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
00751                         do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
00752                         do_fio(&c__1, (char *)&k, (ftnlen)sizeof(integer));
00753                         do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
00754                                 integer));
00755                         do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
00756                                 ;
00757                         do_fio(&c__1, (char *)&jr, (ftnlen)sizeof(integer));
00758                         do_fio(&c__1, (char *)&result[jr], (ftnlen)sizeof(
00759                                 doublereal));
00760                         e_wsfe();
00761                     }
00762 /* L130: */
00763                 }
00764 
00765 L140:
00766                 ;
00767             }
00768 L150:
00769             ;
00770         }
00771 /* L160: */
00772     }
00773 
00774 /*     Summary */
00775 
00776     dlasum_("ZBB", nounit, &nerrs, &ntestt);
00777     return 0;
00778 
00779 
00780 /*     End of ZCHKBB */
00781 
00782 } /* zchkbb_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:17