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 cspt01_(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 complex q__1;
00028
00029
00030 integer i__, j, jc;
00031 real eps;
00032 integer info;
00033 extern logical lsame_(char *, char *);
00034 real anorm;
00035 extern doublereal slamch_(char *);
00036 extern int claset_(char *, integer *, integer *, complex
00037 *, complex *, complex *, integer *);
00038 extern doublereal clansp_(char *, char *, integer *, complex *, real *);
00039 extern int clavsp_(char *, char *, char *, integer *,
00040 integer *, complex *, integer *, complex *, integer *, integer *);
00041 extern doublereal clansy_(char *, char *, integer *, complex *, integer *,
00042 real *);
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 --a;
00117 --afac;
00118 --ipiv;
00119 c_dim1 = *ldc;
00120 c_offset = 1 + c_dim1;
00121 c__ -= c_offset;
00122 --rwork;
00123
00124
00125 if (*n <= 0) {
00126 *resid = 0.f;
00127 return 0;
00128 }
00129
00130
00131
00132 eps = slamch_("Epsilon");
00133 anorm = clansp_("1", uplo, n, &a[1], &rwork[1]);
00134
00135
00136
00137 claset_("Full", n, n, &c_b1, &c_b2, &c__[c_offset], ldc);
00138
00139
00140
00141 clavsp_(uplo, "Transpose", "Non-unit", n, n, &afac[1], &ipiv[1], &c__[
00142 c_offset], ldc, &info);
00143
00144
00145
00146 clavsp_(uplo, "No transpose", "Unit", n, n, &afac[1], &ipiv[1], &c__[
00147 c_offset], ldc, &info);
00148
00149
00150
00151 if (lsame_(uplo, "U")) {
00152 jc = 0;
00153 i__1 = *n;
00154 for (j = 1; j <= i__1; ++j) {
00155 i__2 = j;
00156 for (i__ = 1; i__ <= i__2; ++i__) {
00157 i__3 = i__ + j * c_dim1;
00158 i__4 = i__ + j * c_dim1;
00159 i__5 = jc + i__;
00160 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00161 i__5].i;
00162 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00163
00164 }
00165 jc += j;
00166
00167 }
00168 } else {
00169 jc = 1;
00170 i__1 = *n;
00171 for (j = 1; j <= i__1; ++j) {
00172 i__2 = *n;
00173 for (i__ = j; i__ <= i__2; ++i__) {
00174 i__3 = i__ + j * c_dim1;
00175 i__4 = i__ + j * c_dim1;
00176 i__5 = jc + i__ - j;
00177 q__1.r = c__[i__4].r - a[i__5].r, q__1.i = c__[i__4].i - a[
00178 i__5].i;
00179 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
00180
00181 }
00182 jc = jc + *n - j + 1;
00183
00184 }
00185 }
00186
00187
00188
00189 *resid = clansy_("1", uplo, n, &c__[c_offset], ldc, &rwork[1]);
00190
00191 if (anorm <= 0.f) {
00192 if (*resid != 0.f) {
00193 *resid = 1.f / eps;
00194 }
00195 } else {
00196 *resid = *resid / (real) (*n) / anorm / eps;
00197 }
00198
00199 return 0;
00200
00201
00202
00203 }