zpot03.c
Go to the documentation of this file.
00001 /* zpot03.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 
00020 /* Subroutine */ int zpot03_(char *uplo, integer *n, doublecomplex *a, 
00021         integer *lda, doublecomplex *ainv, integer *ldainv, doublecomplex *
00022         work, integer *ldwork, doublereal *rwork, doublereal *rcond, 
00023         doublereal *resid)
00024 {
00025     /* System generated locals */
00026     integer a_dim1, a_offset, ainv_dim1, ainv_offset, work_dim1, work_offset, 
00027             i__1, i__2, i__3;
00028     doublecomplex z__1;
00029 
00030     /* Builtin functions */
00031     void d_cnjg(doublecomplex *, doublecomplex *);
00032 
00033     /* Local variables */
00034     integer i__, j;
00035     doublereal eps;
00036     extern logical lsame_(char *, char *);
00037     doublereal anorm;
00038     extern /* Subroutine */ int zhemm_(char *, char *, integer *, integer *, 
00039             doublecomplex *, doublecomplex *, integer *, doublecomplex *, 
00040             integer *, doublecomplex *, doublecomplex *, integer *);
00041     extern doublereal dlamch_(char *), zlange_(char *, integer *, 
00042             integer *, doublecomplex *, integer *, doublereal *), 
00043             zlanhe_(char *, char *, integer *, doublecomplex *, integer *, 
00044             doublereal *);
00045     doublereal ainvnm;
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 /*  ZPOT03 computes the residual for a Hermitian matrix times its */
00061 /*  inverse: */
00062 /*     norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ), */
00063 /*  where EPS is the machine epsilon. */
00064 
00065 /*  Arguments */
00066 /*  ========== */
00067 
00068 /*  UPLO    (input) CHARACTER*1 */
00069 /*          Specifies whether the upper or lower triangular part of the */
00070 /*          Hermitian matrix A is stored: */
00071 /*          = 'U':  Upper triangular */
00072 /*          = 'L':  Lower triangular */
00073 
00074 /*  N       (input) INTEGER */
00075 /*          The number of rows and columns of the matrix A.  N >= 0. */
00076 
00077 /*  A       (input) COMPLEX*16 array, dimension (LDA,N) */
00078 /*          The original Hermitian matrix A. */
00079 
00080 /*  LDA     (input) INTEGER */
00081 /*          The leading dimension of the array A.  LDA >= max(1,N) */
00082 
00083 /*  AINV    (input/output) COMPLEX*16 array, dimension (LDAINV,N) */
00084 /*          On entry, the inverse of the matrix A, stored as a Hermitian */
00085 /*          matrix in the same format as A. */
00086 /*          In this version, AINV is expanded into a full matrix and */
00087 /*          multiplied by A, so the opposing triangle of AINV will be */
00088 /*          changed; i.e., if the upper triangular part of AINV is */
00089 /*          stored, the lower triangular part will be used as work space. */
00090 
00091 /*  LDAINV  (input) INTEGER */
00092 /*          The leading dimension of the array AINV.  LDAINV >= max(1,N). */
00093 
00094 /*  WORK    (workspace) COMPLEX*16 array, dimension (LDWORK,N) */
00095 
00096 /*  LDWORK  (input) INTEGER */
00097 /*          The leading dimension of the array WORK.  LDWORK >= max(1,N). */
00098 
00099 /*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N) */
00100 
00101 /*  RCOND   (output) DOUBLE PRECISION */
00102 /*          The reciprocal of the condition number of A, computed as */
00103 /*          ( 1/norm(A) ) / norm(AINV). */
00104 
00105 /*  RESID   (output) DOUBLE PRECISION */
00106 /*          norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS ) */
00107 
00108 /*  ===================================================================== */
00109 
00110 /*     .. Parameters .. */
00111 /*     .. */
00112 /*     .. Local Scalars .. */
00113 /*     .. */
00114 /*     .. External Functions .. */
00115 /*     .. */
00116 /*     .. External Subroutines .. */
00117 /*     .. */
00118 /*     .. Intrinsic Functions .. */
00119 /*     .. */
00120 /*     .. Executable Statements .. */
00121 
00122 /*     Quick exit if N = 0. */
00123 
00124     /* Parameter adjustments */
00125     a_dim1 = *lda;
00126     a_offset = 1 + a_dim1;
00127     a -= a_offset;
00128     ainv_dim1 = *ldainv;
00129     ainv_offset = 1 + ainv_dim1;
00130     ainv -= ainv_offset;
00131     work_dim1 = *ldwork;
00132     work_offset = 1 + work_dim1;
00133     work -= work_offset;
00134     --rwork;
00135 
00136     /* Function Body */
00137     if (*n <= 0) {
00138         *rcond = 1.;
00139         *resid = 0.;
00140         return 0;
00141     }
00142 
00143 /*     Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0. */
00144 
00145     eps = dlamch_("Epsilon");
00146     anorm = zlanhe_("1", uplo, n, &a[a_offset], lda, &rwork[1]);
00147     ainvnm = zlanhe_("1", uplo, n, &ainv[ainv_offset], ldainv, &rwork[1]);
00148     if (anorm <= 0. || ainvnm <= 0.) {
00149         *rcond = 0.;
00150         *resid = 1. / eps;
00151         return 0;
00152     }
00153     *rcond = 1. / anorm / ainvnm;
00154 
00155 /*     Expand AINV into a full matrix and call ZHEMM to multiply */
00156 /*     AINV on the left by A. */
00157 
00158     if (lsame_(uplo, "U")) {
00159         i__1 = *n;
00160         for (j = 1; j <= i__1; ++j) {
00161             i__2 = j - 1;
00162             for (i__ = 1; i__ <= i__2; ++i__) {
00163                 i__3 = j + i__ * ainv_dim1;
00164                 d_cnjg(&z__1, &ainv[i__ + j * ainv_dim1]);
00165                 ainv[i__3].r = z__1.r, ainv[i__3].i = z__1.i;
00166 /* L10: */
00167             }
00168 /* L20: */
00169         }
00170     } else {
00171         i__1 = *n;
00172         for (j = 1; j <= i__1; ++j) {
00173             i__2 = *n;
00174             for (i__ = j + 1; i__ <= i__2; ++i__) {
00175                 i__3 = j + i__ * ainv_dim1;
00176                 d_cnjg(&z__1, &ainv[i__ + j * ainv_dim1]);
00177                 ainv[i__3].r = z__1.r, ainv[i__3].i = z__1.i;
00178 /* L30: */
00179             }
00180 /* L40: */
00181         }
00182     }
00183     z__1.r = -1., z__1.i = -0.;
00184     zhemm_("Left", uplo, n, n, &z__1, &a[a_offset], lda, &ainv[ainv_offset], 
00185             ldainv, &c_b1, &work[work_offset], ldwork);
00186 
00187 /*     Add the identity matrix to WORK . */
00188 
00189     i__1 = *n;
00190     for (i__ = 1; i__ <= i__1; ++i__) {
00191         i__2 = i__ + i__ * work_dim1;
00192         i__3 = i__ + i__ * work_dim1;
00193         z__1.r = work[i__3].r + 1., z__1.i = work[i__3].i + 0.;
00194         work[i__2].r = z__1.r, work[i__2].i = z__1.i;
00195 /* L50: */
00196     }
00197 
00198 /*     Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS) */
00199 
00200     *resid = zlange_("1", n, n, &work[work_offset], ldwork, &rwork[1]);
00201 
00202     *resid = *resid * *rcond / eps / (doublereal) (*n);
00203 
00204     return 0;
00205 
00206 /*     End of ZPOT03 */
00207 
00208 } /* zpot03_ */


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