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
00020 int spocon_(char *uplo, integer *n, real *a, integer *lda,
00021 real *anorm, real *rcond, real *work, integer *iwork, integer *info)
00022 {
00023
00024 integer a_dim1, a_offset, i__1;
00025 real r__1;
00026
00027
00028 integer ix, kase;
00029 real scale;
00030 extern logical lsame_(char *, char *);
00031 integer isave[3];
00032 extern int srscl_(integer *, real *, real *, integer *);
00033 logical upper;
00034 extern int slacn2_(integer *, real *, real *, integer *,
00035 real *, integer *, integer *);
00036 real scalel;
00037 extern doublereal slamch_(char *);
00038 real scaleu;
00039 extern int xerbla_(char *, integer *);
00040 extern integer isamax_(integer *, real *, integer *);
00041 real ainvnm;
00042 char normin[1];
00043 extern int slatrs_(char *, char *, char *, char *,
00044 integer *, real *, integer *, real *, real *, real *, integer *);
00045 real smlnum;
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 a_dim1 = *lda;
00122 a_offset = 1 + a_dim1;
00123 a -= a_offset;
00124 --work;
00125 --iwork;
00126
00127
00128 *info = 0;
00129 upper = lsame_(uplo, "U");
00130 if (! upper && ! lsame_(uplo, "L")) {
00131 *info = -1;
00132 } else if (*n < 0) {
00133 *info = -2;
00134 } else if (*lda < max(1,*n)) {
00135 *info = -4;
00136 } else if (*anorm < 0.f) {
00137 *info = -5;
00138 }
00139 if (*info != 0) {
00140 i__1 = -(*info);
00141 xerbla_("SPOCON", &i__1);
00142 return 0;
00143 }
00144
00145
00146
00147 *rcond = 0.f;
00148 if (*n == 0) {
00149 *rcond = 1.f;
00150 return 0;
00151 } else if (*anorm == 0.f) {
00152 return 0;
00153 }
00154
00155 smlnum = slamch_("Safe minimum");
00156
00157
00158
00159 kase = 0;
00160 *(unsigned char *)normin = 'N';
00161 L10:
00162 slacn2_(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase, isave);
00163 if (kase != 0) {
00164 if (upper) {
00165
00166
00167
00168 slatrs_("Upper", "Transpose", "Non-unit", normin, n, &a[a_offset],
00169 lda, &work[1], &scalel, &work[(*n << 1) + 1], info);
00170 *(unsigned char *)normin = 'Y';
00171
00172
00173
00174 slatrs_("Upper", "No transpose", "Non-unit", normin, n, &a[
00175 a_offset], lda, &work[1], &scaleu, &work[(*n << 1) + 1],
00176 info);
00177 } else {
00178
00179
00180
00181 slatrs_("Lower", "No transpose", "Non-unit", normin, n, &a[
00182 a_offset], lda, &work[1], &scalel, &work[(*n << 1) + 1],
00183 info);
00184 *(unsigned char *)normin = 'Y';
00185
00186
00187
00188 slatrs_("Lower", "Transpose", "Non-unit", normin, n, &a[a_offset],
00189 lda, &work[1], &scaleu, &work[(*n << 1) + 1], info);
00190 }
00191
00192
00193
00194 scale = scalel * scaleu;
00195 if (scale != 1.f) {
00196 ix = isamax_(n, &work[1], &c__1);
00197 if (scale < (r__1 = work[ix], dabs(r__1)) * smlnum || scale ==
00198 0.f) {
00199 goto L20;
00200 }
00201 srscl_(n, &scale, &work[1], &c__1);
00202 }
00203 goto L10;
00204 }
00205
00206
00207
00208 if (ainvnm != 0.f) {
00209 *rcond = 1.f / ainvnm / *anorm;
00210 }
00211
00212 L20:
00213 return 0;
00214
00215
00216
00217 }