zbdt01.c
Go to the documentation of this file.
00001 /* zbdt01.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 doublecomplex c_b7 = {-1.,-0.};
00020 static doublecomplex c_b10 = {1.,0.};
00021 
00022 /* Subroutine */ int zbdt01_(integer *m, integer *n, integer *kd, 
00023         doublecomplex *a, integer *lda, doublecomplex *q, integer *ldq, 
00024         doublereal *d__, doublereal *e, doublecomplex *pt, integer *ldpt, 
00025         doublecomplex *work, doublereal *rwork, doublereal *resid)
00026 {
00027     /* System generated locals */
00028     integer a_dim1, a_offset, pt_dim1, pt_offset, q_dim1, q_offset, i__1, 
00029             i__2, i__3, i__4, i__5, i__6, i__7;
00030     doublereal d__1, d__2;
00031     doublecomplex z__1, z__2, z__3;
00032 
00033     /* Local variables */
00034     integer i__, j;
00035     doublereal eps, anorm;
00036     extern /* Subroutine */ int zgemv_(char *, integer *, integer *, 
00037             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00038             integer *, doublecomplex *, doublecomplex *, integer *), 
00039             zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, 
00040             integer *);
00041     extern doublereal dlamch_(char *), zlange_(char *, integer *, 
00042             integer *, doublecomplex *, integer *, doublereal *), 
00043             dzasum_(integer *, doublecomplex *, integer *);
00044 
00045 
00046 /*  -- LAPACK test routine (version 3.1) -- */
00047 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00048 /*     November 2006 */
00049 
00050 /*     .. Scalar Arguments .. */
00051 /*     .. */
00052 /*     .. Array Arguments .. */
00053 /*     .. */
00054 
00055 /*  Purpose */
00056 /*  ======= */
00057 
00058 /*  ZBDT01 reconstructs a general matrix A from its bidiagonal form */
00059 /*     A = Q * B * P' */
00060 /*  where Q (m by min(m,n)) and P' (min(m,n) by n) are unitary */
00061 /*  matrices and B is bidiagonal. */
00062 
00063 /*  The test ratio to test the reduction is */
00064 /*     RESID = norm( A - Q * B * PT ) / ( n * norm(A) * EPS ) */
00065 /*  where PT = P' and EPS is the machine precision. */
00066 
00067 /*  Arguments */
00068 /*  ========= */
00069 
00070 /*  M       (input) INTEGER */
00071 /*          The number of rows of the matrices A and Q. */
00072 
00073 /*  N       (input) INTEGER */
00074 /*          The number of columns of the matrices A and P'. */
00075 
00076 /*  KD      (input) INTEGER */
00077 /*          If KD = 0, B is diagonal and the array E is not referenced. */
00078 /*          If KD = 1, the reduction was performed by xGEBRD; B is upper */
00079 /*          bidiagonal if M >= N, and lower bidiagonal if M < N. */
00080 /*          If KD = -1, the reduction was performed by xGBBRD; B is */
00081 /*          always upper bidiagonal. */
00082 
00083 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00084 /*          The m by n matrix A. */
00085 
00086 /*  LDA     (input) INTEGER */
00087 /*          The leading dimension of the array A.  LDA >= max(1,M). */
00088 
00089 /*  Q       (input) COMPLEX*16 array, dimension (LDQ,N) */
00090 /*          The m by min(m,n) unitary matrix Q in the reduction */
00091 /*          A = Q * B * P'. */
00092 
00093 /*  LDQ     (input) INTEGER */
00094 /*          The leading dimension of the array Q.  LDQ >= max(1,M). */
00095 
00096 /*  D       (input) DOUBLE PRECISION array, dimension (min(M,N)) */
00097 /*          The diagonal elements of the bidiagonal matrix B. */
00098 
00099 /*  E       (input) DOUBLE PRECISION array, dimension (min(M,N)-1) */
00100 /*          The superdiagonal elements of the bidiagonal matrix B if */
00101 /*          m >= n, or the subdiagonal elements of B if m < n. */
00102 
00103 /*  PT      (input) COMPLEX*16 array, dimension (LDPT,N) */
00104 /*          The min(m,n) by n unitary matrix P' in the reduction */
00105 /*          A = Q * B * P'. */
00106 
00107 /*  LDPT    (input) INTEGER */
00108 /*          The leading dimension of the array PT. */
00109 /*          LDPT >= max(1,min(M,N)). */
00110 
00111 /*  WORK    (workspace) COMPLEX*16 array, dimension (M+N) */
00112 
00113 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (M) */
00114 
00115 /*  RESID   (output) DOUBLE PRECISION */
00116 /*          The test ratio:  norm(A - Q * B * P') / ( n * norm(A) * EPS ) */
00117 
00118 /*  ===================================================================== */
00119 
00120 /*     .. Parameters .. */
00121 /*     .. */
00122 /*     .. Local Scalars .. */
00123 /*     .. */
00124 /*     .. External Functions .. */
00125 /*     .. */
00126 /*     .. External Subroutines .. */
00127 /*     .. */
00128 /*     .. Intrinsic Functions .. */
00129 /*     .. */
00130 /*     .. Executable Statements .. */
00131 
00132 /*     Quick return if possible */
00133 
00134     /* Parameter adjustments */
00135     a_dim1 = *lda;
00136     a_offset = 1 + a_dim1;
00137     a -= a_offset;
00138     q_dim1 = *ldq;
00139     q_offset = 1 + q_dim1;
00140     q -= q_offset;
00141     --d__;
00142     --e;
00143     pt_dim1 = *ldpt;
00144     pt_offset = 1 + pt_dim1;
00145     pt -= pt_offset;
00146     --work;
00147     --rwork;
00148 
00149     /* Function Body */
00150     if (*m <= 0 || *n <= 0) {
00151         *resid = 0.;
00152         return 0;
00153     }
00154 
00155 /*     Compute A - Q * B * P' one column at a time. */
00156 
00157     *resid = 0.;
00158     if (*kd != 0) {
00159 
00160 /*        B is bidiagonal. */
00161 
00162         if (*kd != 0 && *m >= *n) {
00163 
00164 /*           B is upper bidiagonal and M >= N. */
00165 
00166             i__1 = *n;
00167             for (j = 1; j <= i__1; ++j) {
00168                 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00169                 i__2 = *n - 1;
00170                 for (i__ = 1; i__ <= i__2; ++i__) {
00171                     i__3 = *m + i__;
00172                     i__4 = i__;
00173                     i__5 = i__ + j * pt_dim1;
00174                     z__2.r = d__[i__4] * pt[i__5].r, z__2.i = d__[i__4] * pt[
00175                             i__5].i;
00176                     i__6 = i__;
00177                     i__7 = i__ + 1 + j * pt_dim1;
00178                     z__3.r = e[i__6] * pt[i__7].r, z__3.i = e[i__6] * pt[i__7]
00179                             .i;
00180                     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00181                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00182 /* L10: */
00183                 }
00184                 i__2 = *m + *n;
00185                 i__3 = *n;
00186                 i__4 = *n + j * pt_dim1;
00187                 z__1.r = d__[i__3] * pt[i__4].r, z__1.i = d__[i__3] * pt[i__4]
00188                         .i;
00189                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00190                 zgemv_("No transpose", m, n, &c_b7, &q[q_offset], ldq, &work[*
00191                         m + 1], &c__1, &c_b10, &work[1], &c__1);
00192 /* Computing MAX */
00193                 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00194                 *resid = max(d__1,d__2);
00195 /* L20: */
00196             }
00197         } else if (*kd < 0) {
00198 
00199 /*           B is upper bidiagonal and M < N. */
00200 
00201             i__1 = *n;
00202             for (j = 1; j <= i__1; ++j) {
00203                 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00204                 i__2 = *m - 1;
00205                 for (i__ = 1; i__ <= i__2; ++i__) {
00206                     i__3 = *m + i__;
00207                     i__4 = i__;
00208                     i__5 = i__ + j * pt_dim1;
00209                     z__2.r = d__[i__4] * pt[i__5].r, z__2.i = d__[i__4] * pt[
00210                             i__5].i;
00211                     i__6 = i__;
00212                     i__7 = i__ + 1 + j * pt_dim1;
00213                     z__3.r = e[i__6] * pt[i__7].r, z__3.i = e[i__6] * pt[i__7]
00214                             .i;
00215                     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00216                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00217 /* L30: */
00218                 }
00219                 i__2 = *m + *m;
00220                 i__3 = *m;
00221                 i__4 = *m + j * pt_dim1;
00222                 z__1.r = d__[i__3] * pt[i__4].r, z__1.i = d__[i__3] * pt[i__4]
00223                         .i;
00224                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00225                 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00226                         m + 1], &c__1, &c_b10, &work[1], &c__1);
00227 /* Computing MAX */
00228                 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00229                 *resid = max(d__1,d__2);
00230 /* L40: */
00231             }
00232         } else {
00233 
00234 /*           B is lower bidiagonal. */
00235 
00236             i__1 = *n;
00237             for (j = 1; j <= i__1; ++j) {
00238                 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00239                 i__2 = *m + 1;
00240                 i__3 = j * pt_dim1 + 1;
00241                 z__1.r = d__[1] * pt[i__3].r, z__1.i = d__[1] * pt[i__3].i;
00242                 work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00243                 i__2 = *m;
00244                 for (i__ = 2; i__ <= i__2; ++i__) {
00245                     i__3 = *m + i__;
00246                     i__4 = i__ - 1;
00247                     i__5 = i__ - 1 + j * pt_dim1;
00248                     z__2.r = e[i__4] * pt[i__5].r, z__2.i = e[i__4] * pt[i__5]
00249                             .i;
00250                     i__6 = i__;
00251                     i__7 = i__ + j * pt_dim1;
00252                     z__3.r = d__[i__6] * pt[i__7].r, z__3.i = d__[i__6] * pt[
00253                             i__7].i;
00254                     z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i;
00255                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00256 /* L50: */
00257                 }
00258                 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00259                         m + 1], &c__1, &c_b10, &work[1], &c__1);
00260 /* Computing MAX */
00261                 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00262                 *resid = max(d__1,d__2);
00263 /* L60: */
00264             }
00265         }
00266     } else {
00267 
00268 /*        B is diagonal. */
00269 
00270         if (*m >= *n) {
00271             i__1 = *n;
00272             for (j = 1; j <= i__1; ++j) {
00273                 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00274                 i__2 = *n;
00275                 for (i__ = 1; i__ <= i__2; ++i__) {
00276                     i__3 = *m + i__;
00277                     i__4 = i__;
00278                     i__5 = i__ + j * pt_dim1;
00279                     z__1.r = d__[i__4] * pt[i__5].r, z__1.i = d__[i__4] * pt[
00280                             i__5].i;
00281                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00282 /* L70: */
00283                 }
00284                 zgemv_("No transpose", m, n, &c_b7, &q[q_offset], ldq, &work[*
00285                         m + 1], &c__1, &c_b10, &work[1], &c__1);
00286 /* Computing MAX */
00287                 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00288                 *resid = max(d__1,d__2);
00289 /* L80: */
00290             }
00291         } else {
00292             i__1 = *n;
00293             for (j = 1; j <= i__1; ++j) {
00294                 zcopy_(m, &a[j * a_dim1 + 1], &c__1, &work[1], &c__1);
00295                 i__2 = *m;
00296                 for (i__ = 1; i__ <= i__2; ++i__) {
00297                     i__3 = *m + i__;
00298                     i__4 = i__;
00299                     i__5 = i__ + j * pt_dim1;
00300                     z__1.r = d__[i__4] * pt[i__5].r, z__1.i = d__[i__4] * pt[
00301                             i__5].i;
00302                     work[i__3].r = z__1.r, work[i__3].i = z__1.i;
00303 /* L90: */
00304                 }
00305                 zgemv_("No transpose", m, m, &c_b7, &q[q_offset], ldq, &work[*
00306                         m + 1], &c__1, &c_b10, &work[1], &c__1);
00307 /* Computing MAX */
00308                 d__1 = *resid, d__2 = dzasum_(m, &work[1], &c__1);
00309                 *resid = max(d__1,d__2);
00310 /* L100: */
00311             }
00312         }
00313     }
00314 
00315 /*     Compute norm(A - Q * B * P') / ( n * norm(A) * EPS ) */
00316 
00317     anorm = zlange_("1", m, n, &a[a_offset], lda, &rwork[1]);
00318     eps = dlamch_("Precision");
00319 
00320     if (anorm <= 0.) {
00321         if (*resid != 0.) {
00322             *resid = 1. / eps;
00323         }
00324     } else {
00325         if (anorm >= *resid) {
00326             *resid = *resid / anorm / ((doublereal) (*n) * eps);
00327         } else {
00328             if (anorm < 1.) {
00329 /* Computing MIN */
00330                 d__1 = *resid, d__2 = (doublereal) (*n) * anorm;
00331                 *resid = min(d__1,d__2) / anorm / ((doublereal) (*n) * eps);
00332             } else {
00333 /* Computing MIN */
00334                 d__1 = *resid / anorm, d__2 = (doublereal) (*n);
00335                 *resid = min(d__1,d__2) / ((doublereal) (*n) * eps);
00336             }
00337         }
00338     }
00339 
00340     return 0;
00341 
00342 /*     End of ZBDT01 */
00343 
00344 } /* zbdt01_ */


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