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 static complex c_b2 = {1.f,0.f};
00020
00021 int chpt01_(char *uplo, integer *n, complex *a, complex *
00022 afac, integer *ipiv, complex *c__, integer *ldc, real *rwork, real *
00023 resid)
00024 {
00025
00026 integer c_dim1, c_offset, i__1, i__2, i__3, i__4, i__5;
00027 real r__1;
00028 complex q__1;
00029
00030
00031 double r_imag(complex *);
00032
00033
00034 integer i__, j, jc;
00035 real eps;
00036 integer info;
00037 extern logical lsame_(char *, char *);
00038 real anorm;
00039 extern doublereal clanhe_(char *, char *, integer *, complex *, integer *,
00040 real *), clanhp_(char *, char *, integer *,
00041 complex *, real *), slamch_(char *);
00042 extern int claset_(char *, integer *, integer *, complex
00043 *, complex *, complex *, integer *), clavhp_(char *, char
00044 *, char *, integer *, integer *, complex *, integer *, complex *,
00045 integer *, integer *);
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 --a;
00121 --afac;
00122 --ipiv;
00123 c_dim1 = *ldc;
00124 c_offset = 1 + c_dim1;
00125 c__ -= c_offset;
00126 --rwork;
00127
00128
00129 if (*n <= 0) {
00130 *resid = 0.f;
00131 return 0;
00132 }
00133
00134
00135
00136 eps = slamch_("Epsilon");
00137 anorm = clanhp_("1", uplo, n, &a[1], &rwork[1]);
00138
00139
00140
00141
00142 jc = 1;
00143 if (lsame_(uplo, "U")) {
00144 i__1 = *n;
00145 for (j = 1; j <= i__1; ++j) {
00146 if (r_imag(&afac[jc]) != 0.f) {
00147 *resid = 1.f / eps;
00148 return 0;
00149 }
00150 jc = jc + j + 1;
00151
00152 }
00153 } else {
00154 i__1 = *n;
00155 for (j = 1; j <= i__1; ++j) {
00156 if (r_imag(&afac[jc]) != 0.f) {
00157 *resid = 1.f / eps;
00158 return 0;
00159 }
00160 jc = jc + *n - j + 1;
00161
00162 }
00163 }
00164
00165
00166
00167 claset_("Full", n, n, &c_b1, &c_b2, &c__[c_offset], ldc);
00168
00169
00170
00171 clavhp_(uplo, "Conjugate", "Non-unit", n, n, &afac[1], &ipiv[1], &c__[
00172 c_offset], ldc, &info);
00173
00174
00175
00176 clavhp_(uplo, "No transpose", "Unit", n, n, &afac[1], &ipiv[1], &c__[
00177 c_offset], ldc, &info);
00178
00179
00180
00181 if (lsame_(uplo, "U")) {
00182 jc = 0;
00183 i__1 = *n;
00184 for (j = 1; j <= i__1; ++j) {
00185 i__2 = j - 1;
00186 for (i__ = 1; i__ <= i__2; ++i__) {
00187 i__3 = i__ + j * c_dim1;
00188 i__4 = i__ + j * c_dim1;
00189 i__5 = jc + i__;
00190 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00191 i__5].i;
00192 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00193
00194 }
00195 i__2 = j + j * c_dim1;
00196 i__3 = j + j * c_dim1;
00197 i__4 = jc + j;
00198 r__1 = a[i__4].r;
00199 q__1.r = c__[i__3].r - r__1, q__1.i = c__[i__3].i;
00200 c__[i__2].r = q__1.r, c__[i__2].i = q__1.i;
00201 jc += j;
00202
00203 }
00204 } else {
00205 jc = 1;
00206 i__1 = *n;
00207 for (j = 1; j <= i__1; ++j) {
00208 i__2 = j + j * c_dim1;
00209 i__3 = j + j * c_dim1;
00210 i__4 = jc;
00211 r__1 = a[i__4].r;
00212 q__1.r = c__[i__3].r - r__1, q__1.i = c__[i__3].i;
00213 c__[i__2].r = q__1.r, c__[i__2].i = q__1.i;
00214 i__2 = *n;
00215 for (i__ = j + 1; i__ <= i__2; ++i__) {
00216 i__3 = i__ + j * c_dim1;
00217 i__4 = i__ + j * c_dim1;
00218 i__5 = jc + i__ - j;
00219 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00220 i__5].i;
00221 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00222
00223 }
00224 jc = jc + *n - j + 1;
00225
00226 }
00227 }
00228
00229
00230
00231 *resid = clanhe_("1", uplo, n, &c__[c_offset], ldc, &rwork[1]);
00232
00233 if (anorm <= 0.f) {
00234 if (*resid != 0.f) {
00235 *resid = 1.f / eps;
00236 }
00237 } else {
00238 *resid = *resid / (real) (*n) / anorm / eps;
00239 }
00240
00241 return 0;
00242
00243
00244
00245 }