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


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