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 integer c__1 = 1;
00019 static doublereal c_b19 = 1.;
00020
00021 int zppt01_(char *uplo, integer *n, doublecomplex *a,
00022 doublecomplex *afac, doublereal *rwork, doublereal *resid)
00023 {
00024
00025 integer i__1, i__2, i__3, i__4, i__5;
00026 doublereal d__1;
00027 doublecomplex z__1;
00028
00029
00030 double d_imag(doublecomplex *);
00031
00032
00033 integer i__, k, kc;
00034 doublecomplex tc;
00035 doublereal tr, eps;
00036 extern int zhpr_(char *, integer *, doublereal *,
00037 doublecomplex *, integer *, doublecomplex *);
00038 extern logical lsame_(char *, char *);
00039 doublereal anorm;
00040 extern int zscal_(integer *, doublecomplex *,
00041 doublecomplex *, integer *);
00042 extern VOID zdotc_(doublecomplex *, integer *,
00043 doublecomplex *, integer *, doublecomplex *, integer *);
00044 extern int ztpmv_(char *, char *, char *, integer *,
00045 doublecomplex *, doublecomplex *, integer *);
00046 extern doublereal dlamch_(char *), zlanhp_(char *, char *,
00047 integer *, doublecomplex *, doublereal *);
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 --rwork;
00115 --afac;
00116 --a;
00117
00118
00119 if (*n <= 0) {
00120 *resid = 0.;
00121 return 0;
00122 }
00123
00124
00125
00126 eps = dlamch_("Epsilon");
00127 anorm = zlanhp_("1", uplo, n, &a[1], &rwork[1]);
00128 if (anorm <= 0.) {
00129 *resid = 1. / eps;
00130 return 0;
00131 }
00132
00133
00134
00135
00136 kc = 1;
00137 if (lsame_(uplo, "U")) {
00138 i__1 = *n;
00139 for (k = 1; k <= i__1; ++k) {
00140 if (d_imag(&afac[kc]) != 0.) {
00141 *resid = 1. / eps;
00142 return 0;
00143 }
00144 kc = kc + k + 1;
00145
00146 }
00147 } else {
00148 i__1 = *n;
00149 for (k = 1; k <= i__1; ++k) {
00150 if (d_imag(&afac[kc]) != 0.) {
00151 *resid = 1. / eps;
00152 return 0;
00153 }
00154 kc = kc + *n - k + 1;
00155
00156 }
00157 }
00158
00159
00160
00161 if (lsame_(uplo, "U")) {
00162 kc = *n * (*n - 1) / 2 + 1;
00163 for (k = *n; k >= 1; --k) {
00164
00165
00166
00167 zdotc_(&z__1, &k, &afac[kc], &c__1, &afac[kc], &c__1);
00168 tr = z__1.r;
00169 i__1 = kc + k - 1;
00170 afac[i__1].r = tr, afac[i__1].i = 0.;
00171
00172
00173
00174 if (k > 1) {
00175 i__1 = k - 1;
00176 ztpmv_("Upper", "Conjugate", "Non-unit", &i__1, &afac[1], &
00177 afac[kc], &c__1);
00178 kc -= k - 1;
00179 }
00180
00181 }
00182
00183
00184
00185 kc = 1;
00186 i__1 = *n;
00187 for (k = 1; k <= i__1; ++k) {
00188 i__2 = k - 1;
00189 for (i__ = 1; i__ <= i__2; ++i__) {
00190 i__3 = kc + i__ - 1;
00191 i__4 = kc + i__ - 1;
00192 i__5 = kc + i__ - 1;
00193 z__1.r = afac[i__4].r - a[i__5].r, z__1.i = afac[i__4].i - a[
00194 i__5].i;
00195 afac[i__3].r = z__1.r, afac[i__3].i = z__1.i;
00196
00197 }
00198 i__2 = kc + k - 1;
00199 i__3 = kc + k - 1;
00200 i__4 = kc + k - 1;
00201 d__1 = a[i__4].r;
00202 z__1.r = afac[i__3].r - d__1, z__1.i = afac[i__3].i;
00203 afac[i__2].r = z__1.r, afac[i__2].i = z__1.i;
00204 kc += k;
00205
00206 }
00207
00208
00209
00210 } else {
00211 kc = *n * (*n + 1) / 2;
00212 for (k = *n; k >= 1; --k) {
00213
00214
00215
00216
00217 if (k < *n) {
00218 i__1 = *n - k;
00219 zhpr_("Lower", &i__1, &c_b19, &afac[kc + 1], &c__1, &afac[kc
00220 + *n - k + 1]);
00221 }
00222
00223
00224
00225 i__1 = kc;
00226 tc.r = afac[i__1].r, tc.i = afac[i__1].i;
00227 i__1 = *n - k + 1;
00228 zscal_(&i__1, &tc, &afac[kc], &c__1);
00229
00230 kc -= *n - k + 2;
00231
00232 }
00233
00234
00235
00236 kc = 1;
00237 i__1 = *n;
00238 for (k = 1; k <= i__1; ++k) {
00239 i__2 = kc;
00240 i__3 = kc;
00241 i__4 = kc;
00242 d__1 = a[i__4].r;
00243 z__1.r = afac[i__3].r - d__1, z__1.i = afac[i__3].i;
00244 afac[i__2].r = z__1.r, afac[i__2].i = z__1.i;
00245 i__2 = *n;
00246 for (i__ = k + 1; i__ <= i__2; ++i__) {
00247 i__3 = kc + i__ - k;
00248 i__4 = kc + i__ - k;
00249 i__5 = kc + i__ - k;
00250 z__1.r = afac[i__4].r - a[i__5].r, z__1.i = afac[i__4].i - a[
00251 i__5].i;
00252 afac[i__3].r = z__1.r, afac[i__3].i = z__1.i;
00253
00254 }
00255 kc = kc + *n - k + 1;
00256
00257 }
00258 }
00259
00260
00261
00262 *resid = zlanhp_("1", uplo, n, &afac[1], &rwork[1]);
00263
00264 *resid = *resid / (doublereal) (*n) / anorm / eps;
00265
00266 return 0;
00267
00268
00269
00270 }