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 = {1.,0.};
00019 static doublecomplex c_b2 = {-1.,0.};
00020 static integer c__1 = 1;
00021
00022 int zgetrf_(integer *m, integer *n, doublecomplex *a,
00023 integer *lda, integer *ipiv, integer *info)
00024 {
00025
00026 integer a_dim1, a_offset, i__1, i__2, i__3;
00027 doublecomplex z__1;
00028
00029
00030 double z_abs(doublecomplex *);
00031 void z_div(doublecomplex *, doublecomplex *, doublecomplex *);
00032
00033
00034 integer i__, j, ipivstart, jpivstart, jp;
00035 doublecomplex tmp;
00036 integer kcols;
00037 doublereal sfmin;
00038 extern int zscal_(integer *, doublecomplex *,
00039 doublecomplex *, integer *), zgemm_(char *, char *, integer *,
00040 integer *, integer *, doublecomplex *, doublecomplex *, integer *,
00041 doublecomplex *, integer *, doublecomplex *, doublecomplex *,
00042 integer *);
00043 integer nstep;
00044 extern int ztrsm_(char *, char *, char *, char *,
00045 integer *, integer *, doublecomplex *, doublecomplex *, integer *,
00046 doublecomplex *, integer *);
00047 integer kahead;
00048 extern doublereal dlamch_(char *);
00049 extern logical disnan_(doublereal *);
00050 extern int xerbla_(char *, integer *);
00051 doublereal pivmag;
00052 integer npived;
00053 extern integer izamax_(integer *, doublecomplex *, integer *);
00054 integer kstart, ntopiv;
00055 extern int zlaswp_(integer *, doublecomplex *, integer *,
00056 integer *, integer *, integer *, integer *);
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
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 a_dim1 = *lda;
00163 a_offset = 1 + a_dim1;
00164 a -= a_offset;
00165 --ipiv;
00166
00167
00168 *info = 0;
00169 if (*m < 0) {
00170 *info = -1;
00171 } else if (*n < 0) {
00172 *info = -2;
00173 } else if (*lda < max(1,*m)) {
00174 *info = -4;
00175 }
00176 if (*info != 0) {
00177 i__1 = -(*info);
00178 xerbla_("ZGETRF", &i__1);
00179 return 0;
00180 }
00181
00182
00183
00184 if (*m == 0 || *n == 0) {
00185 return 0;
00186 }
00187
00188
00189
00190 sfmin = dlamch_("S");
00191
00192 nstep = min(*m,*n);
00193 i__1 = nstep;
00194 for (j = 1; j <= i__1; ++j) {
00195 kahead = j & -j;
00196 kstart = j + 1 - kahead;
00197
00198 i__2 = kahead, i__3 = *m - j;
00199 kcols = min(i__2,i__3);
00200
00201
00202
00203 i__2 = *m - j + 1;
00204 jp = j - 1 + izamax_(&i__2, &a[j + j * a_dim1], &c__1);
00205 ipiv[j] = jp;
00206
00207 if (jp != j) {
00208 i__2 = j + j * a_dim1;
00209 tmp.r = a[i__2].r, tmp.i = a[i__2].i;
00210 i__2 = j + j * a_dim1;
00211 i__3 = jp + j * a_dim1;
00212 a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
00213 i__2 = jp + j * a_dim1;
00214 a[i__2].r = tmp.r, a[i__2].i = tmp.i;
00215 }
00216
00217 ntopiv = 1;
00218 ipivstart = j;
00219 jpivstart = j - ntopiv;
00220 while(ntopiv < kahead) {
00221 zlaswp_(&ntopiv, &a[jpivstart * a_dim1 + 1], lda, &ipivstart, &j,
00222 &ipiv[1], &c__1);
00223 ipivstart -= ntopiv;
00224 ntopiv <<= 1;
00225 jpivstart -= ntopiv;
00226 }
00227
00228 zlaswp_(&kcols, &a[(j + 1) * a_dim1 + 1], lda, &kstart, &j, &ipiv[1],
00229 &c__1);
00230
00231 pivmag = z_abs(&a[j + j * a_dim1]);
00232 if (pivmag != 0. && ! disnan_(&pivmag)) {
00233 if (pivmag >= sfmin) {
00234 i__2 = *m - j;
00235 z_div(&z__1, &c_b1, &a[j + j * a_dim1]);
00236 zscal_(&i__2, &z__1, &a[j + 1 + j * a_dim1], &c__1);
00237 } else {
00238 i__2 = *m - j;
00239 for (i__ = 1; i__ <= i__2; ++i__) {
00240 i__3 = j + i__ + j * a_dim1;
00241 z_div(&z__1, &a[j + i__ + j * a_dim1], &a[j + j * a_dim1])
00242 ;
00243 a[i__3].r = z__1.r, a[i__3].i = z__1.i;
00244 }
00245 }
00246 } else if (pivmag == 0. && *info == 0) {
00247 *info = j;
00248 }
00249
00250 ztrsm_("Left", "Lower", "No transpose", "Unit", &kahead, &kcols, &
00251 c_b1, &a[kstart + kstart * a_dim1], lda, &a[kstart + (j + 1) *
00252 a_dim1], lda);
00253
00254 i__2 = *m - j;
00255 zgemm_("No transpose", "No transpose", &i__2, &kcols, &kahead, &c_b2,
00256 &a[j + 1 + kstart * a_dim1], lda, &a[kstart + (j + 1) *
00257 a_dim1], lda, &c_b1, &a[j + 1 + (j + 1) * a_dim1], lda);
00258 }
00259
00260 npived = nstep & -nstep;
00261 j = nstep - npived;
00262 while(j > 0) {
00263 ntopiv = j & -j;
00264 i__1 = j + 1;
00265 zlaswp_(&ntopiv, &a[(j - ntopiv + 1) * a_dim1 + 1], lda, &i__1, &
00266 nstep, &ipiv[1], &c__1);
00267 j -= ntopiv;
00268 }
00269
00270 if (*m < *n) {
00271 i__1 = *n - *m;
00272 zlaswp_(&i__1, &a[(*m + kcols + 1) * a_dim1 + 1], lda, &c__1, m, &
00273 ipiv[1], &c__1);
00274 i__1 = *n - *m;
00275 ztrsm_("Left", "Lower", "No transpose", "Unit", m, &i__1, &c_b1, &a[
00276 a_offset], lda, &a[(*m + kcols + 1) * a_dim1 + 1], lda);
00277 }
00278 return 0;
00279
00280
00281
00282 }