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 doublecomplex c_b1 = {0.,0.};
00019
00020 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
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
00031 void d_cnjg(doublecomplex *, doublecomplex *);
00032
00033
00034 integer i__, j;
00035 doublereal eps;
00036 extern logical lsame_(char *, char *);
00037 doublereal anorm;
00038 extern 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
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
00123
00124
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
00137 if (*n <= 0) {
00138 *rcond = 1.;
00139 *resid = 0.;
00140 return 0;
00141 }
00142
00143
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
00156
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
00167 }
00168
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
00179 }
00180
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
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
00196 }
00197
00198
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
00207
00208 }