zpbt01.c
Go to the documentation of this file.
00001 /* zpbt01.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 doublereal c_b17 = 1.;
00020 
00021 /* Subroutine */ int zpbt01_(char *uplo, integer *n, integer *kd, 
00022         doublecomplex *a, integer *lda, doublecomplex *afac, integer *ldafac, 
00023         doublereal *rwork, doublereal *resid)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, afac_dim1, afac_offset, i__1, i__2, i__3, i__4, 
00027             i__5;
00028     doublecomplex z__1;
00029 
00030     /* Builtin functions */
00031     double d_imag(doublecomplex *);
00032 
00033     /* Local variables */
00034     integer i__, j, k, kc, ml, mu;
00035     doublereal akk, eps;
00036     integer klen;
00037     extern /* Subroutine */ int zher_(char *, integer *, doublereal *, 
00038             doublecomplex *, integer *, doublecomplex *, integer *);
00039     extern logical lsame_(char *, char *);
00040     doublereal anorm;
00041     extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, 
00042             doublecomplex *, integer *, doublecomplex *, integer *);
00043     extern /* Subroutine */ int ztrmv_(char *, char *, char *, integer *, 
00044             doublecomplex *, integer *, doublecomplex *, integer *);
00045     extern doublereal dlamch_(char *), zlanhb_(char *, char *, 
00046             integer *, integer *, doublecomplex *, integer *, doublereal *);
00047     extern /* Subroutine */ int zdscal_(integer *, doublereal *, 
00048             doublecomplex *, integer *);
00049 
00050 
00051 /*  -- LAPACK test routine (version 3.1) -- */
00052 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00053 /*     November 2006 */
00054 
00055 /*     .. Scalar Arguments .. */
00056 /*     .. */
00057 /*     .. Array Arguments .. */
00058 /*     .. */
00059 
00060 /*  Purpose */
00061 /*  ======= */
00062 
00063 /*  ZPBT01 reconstructs a Hermitian positive definite band matrix A from */
00064 /*  its L*L' or U'*U factorization and computes the residual */
00065 /*     norm( L*L' - A ) / ( N * norm(A) * EPS ) or */
00066 /*     norm( U'*U - A ) / ( N * norm(A) * EPS ), */
00067 /*  where EPS is the machine epsilon, L' is the conjugate transpose of */
00068 /*  L, and U' is the conjugate transpose of U. */
00069 
00070 /*  Arguments */
00071 /*  ========= */
00072 
00073 /*  UPLO    (input) CHARACTER*1 */
00074 /*          Specifies whether the upper or lower triangular part of the */
00075 /*          Hermitian matrix A is stored: */
00076 /*          = 'U':  Upper triangular */
00077 /*          = 'L':  Lower triangular */
00078 
00079 /*  N       (input) INTEGER */
00080 /*          The number of rows and columns of the matrix A.  N >= 0. */
00081 
00082 /*  KD      (input) INTEGER */
00083 /*          The number of super-diagonals of the matrix A if UPLO = 'U', */
00084 /*          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0. */
00085 
00086 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00087 /*          The original Hermitian band matrix A.  If UPLO = 'U', the */
00088 /*          upper triangular part of A is stored as a band matrix; if */
00089 /*          UPLO = 'L', the lower triangular part of A is stored.  The */
00090 /*          columns of the appropriate triangle are stored in the columns */
00091 /*          of A and the diagonals of the triangle are stored in the rows */
00092 /*          of A.  See ZPBTRF for further details. */
00093 
00094 /*  LDA     (input) INTEGER. */
00095 /*          The leading dimension of the array A.  LDA >= max(1,KD+1). */
00096 
00097 /*  AFAC    (input) COMPLEX*16 array, dimension (LDAFAC,N) */
00098 /*          The factored form of the matrix A.  AFAC contains the factor */
00099 /*          L or U from the L*L' or U'*U factorization in band storage */
00100 /*          format, as computed by ZPBTRF. */
00101 
00102 /*  LDAFAC  (input) INTEGER */
00103 /*          The leading dimension of the array AFAC. */
00104 /*          LDAFAC >= max(1,KD+1). */
00105 
00106 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00107 
00108 /*  RESID   (output) DOUBLE PRECISION */
00109 /*          If UPLO = 'L', norm(L*L' - A) / ( N * norm(A) * EPS ) */
00110 /*          If UPLO = 'U', norm(U'*U - A) / ( N * norm(A) * EPS ) */
00111 
00112 /*  ===================================================================== */
00113 
00114 
00115 /*     .. Parameters .. */
00116 /*     .. */
00117 /*     .. Local Scalars .. */
00118 /*     .. */
00119 /*     .. External Functions .. */
00120 /*     .. */
00121 /*     .. External Subroutines .. */
00122 /*     .. */
00123 /*     .. Intrinsic Functions .. */
00124 /*     .. */
00125 /*     .. Executable Statements .. */
00126 
00127 /*     Quick exit if N = 0. */
00128 
00129     /* Parameter adjustments */
00130     a_dim1 = *lda;
00131     a_offset = 1 + a_dim1;
00132     a -= a_offset;
00133     afac_dim1 = *ldafac;
00134     afac_offset = 1 + afac_dim1;
00135     afac -= afac_offset;
00136     --rwork;
00137 
00138     /* Function Body */
00139     if (*n <= 0) {
00140         *resid = 0.;
00141         return 0;
00142     }
00143 
00144 /*     Exit with RESID = 1/EPS if ANORM = 0. */
00145 
00146     eps = dlamch_("Epsilon");
00147     anorm = zlanhb_("1", uplo, n, kd, &a[a_offset], lda, &rwork[1]);
00148     if (anorm <= 0.) {
00149         *resid = 1. / eps;
00150         return 0;
00151     }
00152 
00153 /*     Check the imaginary parts of the diagonal elements and return with */
00154 /*     an error code if any are nonzero. */
00155 
00156     if (lsame_(uplo, "U")) {
00157         i__1 = *n;
00158         for (j = 1; j <= i__1; ++j) {
00159             if (d_imag(&afac[*kd + 1 + j * afac_dim1]) != 0.) {
00160                 *resid = 1. / eps;
00161                 return 0;
00162             }
00163 /* L10: */
00164         }
00165     } else {
00166         i__1 = *n;
00167         for (j = 1; j <= i__1; ++j) {
00168             if (d_imag(&afac[j * afac_dim1 + 1]) != 0.) {
00169                 *resid = 1. / eps;
00170                 return 0;
00171             }
00172 /* L20: */
00173         }
00174     }
00175 
00176 /*     Compute the product U'*U, overwriting U. */
00177 
00178     if (lsame_(uplo, "U")) {
00179         for (k = *n; k >= 1; --k) {
00180 /* Computing MAX */
00181             i__1 = 1, i__2 = *kd + 2 - k;
00182             kc = max(i__1,i__2);
00183             klen = *kd + 1 - kc;
00184 
00185 /*           Compute the (K,K) element of the result. */
00186 
00187             i__1 = klen + 1;
00188             zdotc_(&z__1, &i__1, &afac[kc + k * afac_dim1], &c__1, &afac[kc + 
00189                     k * afac_dim1], &c__1);
00190             akk = z__1.r;
00191             i__1 = *kd + 1 + k * afac_dim1;
00192             afac[i__1].r = akk, afac[i__1].i = 0.;
00193 
00194 /*           Compute the rest of column K. */
00195 
00196             if (klen > 0) {
00197                 i__1 = *ldafac - 1;
00198                 ztrmv_("Upper", "Conjugate", "Non-unit", &klen, &afac[*kd + 1 
00199                         + (k - klen) * afac_dim1], &i__1, &afac[kc + k * 
00200                         afac_dim1], &c__1);
00201             }
00202 
00203 /* L30: */
00204         }
00205 
00206 /*     UPLO = 'L':  Compute the product L*L', overwriting L. */
00207 
00208     } else {
00209         for (k = *n; k >= 1; --k) {
00210 /* Computing MIN */
00211             i__1 = *kd, i__2 = *n - k;
00212             klen = min(i__1,i__2);
00213 
00214 /*           Add a multiple of column K of the factor L to each of */
00215 /*           columns K+1 through N. */
00216 
00217             if (klen > 0) {
00218                 i__1 = *ldafac - 1;
00219                 zher_("Lower", &klen, &c_b17, &afac[k * afac_dim1 + 2], &c__1, 
00220                          &afac[(k + 1) * afac_dim1 + 1], &i__1);
00221             }
00222 
00223 /*           Scale column K by the diagonal element. */
00224 
00225             i__1 = k * afac_dim1 + 1;
00226             akk = afac[i__1].r;
00227             i__1 = klen + 1;
00228             zdscal_(&i__1, &akk, &afac[k * afac_dim1 + 1], &c__1);
00229 
00230 /* L40: */
00231         }
00232     }
00233 
00234 /*     Compute the difference  L*L' - A  or  U'*U - A. */
00235 
00236     if (lsame_(uplo, "U")) {
00237         i__1 = *n;
00238         for (j = 1; j <= i__1; ++j) {
00239 /* Computing MAX */
00240             i__2 = 1, i__3 = *kd + 2 - j;
00241             mu = max(i__2,i__3);
00242             i__2 = *kd + 1;
00243             for (i__ = mu; i__ <= i__2; ++i__) {
00244                 i__3 = i__ + j * afac_dim1;
00245                 i__4 = i__ + j * afac_dim1;
00246                 i__5 = i__ + j * a_dim1;
00247                 z__1.r = afac[i__4].r - a[i__5].r, z__1.i = afac[i__4].i - a[
00248                         i__5].i;
00249                 afac[i__3].r = z__1.r, afac[i__3].i = z__1.i;
00250 /* L50: */
00251             }
00252 /* L60: */
00253         }
00254     } else {
00255         i__1 = *n;
00256         for (j = 1; j <= i__1; ++j) {
00257 /* Computing MIN */
00258             i__2 = *kd + 1, i__3 = *n - j + 1;
00259             ml = min(i__2,i__3);
00260             i__2 = ml;
00261             for (i__ = 1; i__ <= i__2; ++i__) {
00262                 i__3 = i__ + j * afac_dim1;
00263                 i__4 = i__ + j * afac_dim1;
00264                 i__5 = i__ + j * a_dim1;
00265                 z__1.r = afac[i__4].r - a[i__5].r, z__1.i = afac[i__4].i - a[
00266                         i__5].i;
00267                 afac[i__3].r = z__1.r, afac[i__3].i = z__1.i;
00268 /* L70: */
00269             }
00270 /* L80: */
00271         }
00272     }
00273 
00274 /*     Compute norm( L*L' - A ) / ( N * norm(A) * EPS ) */
00275 
00276     *resid = zlanhb_("1", uplo, n, kd, &afac[afac_offset], ldafac, &rwork[1]);
00277 
00278     *resid = *resid / (doublereal) (*n) / anorm / eps;
00279 
00280     return 0;
00281 
00282 /*     End of ZPBT01 */
00283 
00284 } /* zpbt01_ */


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