chpt01.c
Go to the documentation of this file.
00001 /* chpt01.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 
00021 /* Subroutine */ int chpt01_(char *uplo, integer *n, complex *a, complex *
00022         afac, integer *ipiv, complex *c__, integer *ldc, real *rwork, real *
00023         resid)
00024 {
00025     /* System generated locals */
00026     integer c_dim1, c_offset, i__1, i__2, i__3, i__4, i__5;
00027     real r__1;
00028     complex q__1;
00029 
00030     /* Builtin functions */
00031     double r_imag(complex *);
00032 
00033     /* Local variables */
00034     integer i__, j, jc;
00035     real eps;
00036     integer info;
00037     extern logical lsame_(char *, char *);
00038     real anorm;
00039     extern doublereal clanhe_(char *, char *, integer *, complex *, integer *, 
00040              real *), clanhp_(char *, char *, integer *, 
00041             complex *, real *), slamch_(char *);
00042     extern /* Subroutine */ int claset_(char *, integer *, integer *, complex 
00043             *, complex *, complex *, integer *), clavhp_(char *, char 
00044             *, char *, integer *, integer *, complex *, integer *, complex *, 
00045             integer *, integer *);
00046 
00047 
00048 /*  -- LAPACK test routine (version 3.1) -- */
00049 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00050 /*     November 2006 */
00051 
00052 /*     .. Scalar Arguments .. */
00053 /*     .. */
00054 /*     .. Array Arguments .. */
00055 /*     .. */
00056 
00057 /*  Purpose */
00058 /*  ======= */
00059 
00060 /*  CHPT01 reconstructs a Hermitian indefinite packed matrix A from its */
00061 /*  block L*D*L' or U*D*U' factorization and computes the residual */
00062 /*     norm( C - A ) / ( N * norm(A) * EPS ), */
00063 /*  where C is the reconstructed matrix, EPS is the machine epsilon, */
00064 /*  L' is the conjugate transpose of L, and U' is the conjugate transpose */
00065 /*  of U. */
00066 
00067 /*  Arguments */
00068 /*  ========== */
00069 
00070 /*  UPLO    (input) CHARACTER*1 */
00071 /*          Specifies whether the upper or lower triangular part of the */
00072 /*          Hermitian matrix A is stored: */
00073 /*          = 'U':  Upper triangular */
00074 /*          = 'L':  Lower triangular */
00075 
00076 /*  N       (input) INTEGER */
00077 /*          The number of rows and columns of the matrix A.  N >= 0. */
00078 
00079 /*  A       (input) COMPLEX array, dimension (N*(N+1)/2) */
00080 /*          The original Hermitian matrix A, stored as a packed */
00081 /*          triangular matrix. */
00082 
00083 /*  AFAC    (input) COMPLEX array, dimension (N*(N+1)/2) */
00084 /*          The factored form of the matrix A, stored as a packed */
00085 /*          triangular matrix.  AFAC contains the block diagonal matrix D */
00086 /*          and the multipliers used to obtain the factor L or U from the */
00087 /*          block L*D*L' or U*D*U' factorization as computed by CHPTRF. */
00088 
00089 /*  IPIV    (input) INTEGER array, dimension (N) */
00090 /*          The pivot indices from CHPTRF. */
00091 
00092 /*  C       (workspace) COMPLEX array, dimension (LDC,N) */
00093 
00094 /*  LDC     (integer) INTEGER */
00095 /*          The leading dimension of the array C.  LDC >= max(1,N). */
00096 
00097 /*  RWORK   (workspace) REAL array, dimension (N) */
00098 
00099 /*  RESID   (output) REAL */
00100 /*          If UPLO = 'L', norm(L*D*L' - A) / ( N * norm(A) * EPS ) */
00101 /*          If UPLO = 'U', norm(U*D*U' - A) / ( N * norm(A) * EPS ) */
00102 
00103 /*  ===================================================================== */
00104 
00105 /*     .. Parameters .. */
00106 /*     .. */
00107 /*     .. Local Scalars .. */
00108 /*     .. */
00109 /*     .. External Functions .. */
00110 /*     .. */
00111 /*     .. External Subroutines .. */
00112 /*     .. */
00113 /*     .. Intrinsic Functions .. */
00114 /*     .. */
00115 /*     .. Executable Statements .. */
00116 
00117 /*     Quick exit if N = 0. */
00118 
00119     /* Parameter adjustments */
00120     --a;
00121     --afac;
00122     --ipiv;
00123     c_dim1 = *ldc;
00124     c_offset = 1 + c_dim1;
00125     c__ -= c_offset;
00126     --rwork;
00127 
00128     /* Function Body */
00129     if (*n <= 0) {
00130         *resid = 0.f;
00131         return 0;
00132     }
00133 
00134 /*     Determine EPS and the norm of A. */
00135 
00136     eps = slamch_("Epsilon");
00137     anorm = clanhp_("1", uplo, n, &a[1], &rwork[1]);
00138 
00139 /*     Check the imaginary parts of the diagonal elements and return with */
00140 /*     an error code if any are nonzero. */
00141 
00142     jc = 1;
00143     if (lsame_(uplo, "U")) {
00144         i__1 = *n;
00145         for (j = 1; j <= i__1; ++j) {
00146             if (r_imag(&afac[jc]) != 0.f) {
00147                 *resid = 1.f / eps;
00148                 return 0;
00149             }
00150             jc = jc + j + 1;
00151 /* L10: */
00152         }
00153     } else {
00154         i__1 = *n;
00155         for (j = 1; j <= i__1; ++j) {
00156             if (r_imag(&afac[jc]) != 0.f) {
00157                 *resid = 1.f / eps;
00158                 return 0;
00159             }
00160             jc = jc + *n - j + 1;
00161 /* L20: */
00162         }
00163     }
00164 
00165 /*     Initialize C to the identity matrix. */
00166 
00167     claset_("Full", n, n, &c_b1, &c_b2, &c__[c_offset], ldc);
00168 
00169 /*     Call CLAVHP to form the product D * U' (or D * L' ). */
00170 
00171     clavhp_(uplo, "Conjugate", "Non-unit", n, n, &afac[1], &ipiv[1], &c__[
00172             c_offset], ldc, &info);
00173 
00174 /*     Call CLAVHP again to multiply by U ( or L ). */
00175 
00176     clavhp_(uplo, "No transpose", "Unit", n, n, &afac[1], &ipiv[1], &c__[
00177             c_offset], ldc, &info);
00178 
00179 /*     Compute the difference  C - A . */
00180 
00181     if (lsame_(uplo, "U")) {
00182         jc = 0;
00183         i__1 = *n;
00184         for (j = 1; j <= i__1; ++j) {
00185             i__2 = j - 1;
00186             for (i__ = 1; i__ <= i__2; ++i__) {
00187                 i__3 = i__ + j * c_dim1;
00188                 i__4 = i__ + j * c_dim1;
00189                 i__5 = jc + i__;
00190                 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00191                         i__5].i;
00192                 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00193 /* L30: */
00194             }
00195             i__2 = j + j * c_dim1;
00196             i__3 = j + j * c_dim1;
00197             i__4 = jc + j;
00198             r__1 = a[i__4].r;
00199             q__1.r = c__[i__3].r - r__1, q__1.i = c__[i__3].i;
00200             c__[i__2].r = q__1.r, c__[i__2].i = q__1.i;
00201             jc += j;
00202 /* L40: */
00203         }
00204     } else {
00205         jc = 1;
00206         i__1 = *n;
00207         for (j = 1; j <= i__1; ++j) {
00208             i__2 = j + j * c_dim1;
00209             i__3 = j + j * c_dim1;
00210             i__4 = jc;
00211             r__1 = a[i__4].r;
00212             q__1.r = c__[i__3].r - r__1, q__1.i = c__[i__3].i;
00213             c__[i__2].r = q__1.r, c__[i__2].i = q__1.i;
00214             i__2 = *n;
00215             for (i__ = j + 1; i__ <= i__2; ++i__) {
00216                 i__3 = i__ + j * c_dim1;
00217                 i__4 = i__ + j * c_dim1;
00218                 i__5 = jc + i__ - j;
00219                 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00220                         i__5].i;
00221                 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00222 /* L50: */
00223             }
00224             jc = jc + *n - j + 1;
00225 /* L60: */
00226         }
00227     }
00228 
00229 /*     Compute norm( C - A ) / ( N * norm(A) * EPS ) */
00230 
00231     *resid = clanhe_("1", uplo, n, &c__[c_offset], ldc, &rwork[1]);
00232 
00233     if (anorm <= 0.f) {
00234         if (*resid != 0.f) {
00235             *resid = 1.f / eps;
00236         }
00237     } else {
00238         *resid = *resid / (real) (*n) / anorm / eps;
00239     }
00240 
00241     return 0;
00242 
00243 /*     End of CHPT01 */
00244 
00245 } /* chpt01_ */


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