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


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