00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015
00016
00017
00018 static complex c_b1 = {0.f,0.f};
00019
00020 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
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
00030 integer i__, j;
00031 real eps;
00032 extern logical lsame_(char *, char *);
00033 real anorm;
00034 extern 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
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
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
00134 if (*n <= 0) {
00135 *rcond = 1.f;
00136 *resid = 0.f;
00137 return 0;
00138 }
00139
00140
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
00153
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
00164 }
00165
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
00176 }
00177
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
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
00193 }
00194
00195
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
00204
00205 }