00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "amesos_cholmod_template.h"
00020
00021 #ifdef MASK
00022 static int TEMPLATE (cholmod_rowfac_mask)
00023 #else
00024 static int TEMPLATE (cholmod_rowfac)
00025 #endif
00026 (
00027
00028 cholmod_sparse *A,
00029 cholmod_sparse *F,
00030 double beta [2],
00031 size_t kstart,
00032 size_t kend,
00033 #ifdef MASK
00034
00035 Int *mask,
00036 Int *RLinkUp,
00037 #endif
00038
00039 cholmod_factor *L,
00040
00041 cholmod_common *Common
00042 )
00043 {
00044 double yx [2], lx [2], fx [2], dk [1], di [1], fl = 0 ;
00045 #ifdef ZOMPLEX
00046 double yz [1], lz [1], fz [1] ;
00047 #endif
00048 double *Ax, *Az, *Lx, *Lz, *Wx, *Wz, *Fx, *Fz ;
00049 Int *Ap, *Anz, *Ai, *Lp, *Lnz, *Li, *Lnext, *Flag, *Stack, *Fp, *Fi, *Fnz,
00050 *Iwork ;
00051 Int i, p, k, t, pf, pfend, top, s, mark, pend, n, lnz, is_ll, multadds,
00052 use_dbound, packed, stype, Fpacked, sorted, nzmax, len, parent ;
00053 #ifndef REAL
00054 Int dk_imaginary ;
00055 #endif
00056
00057
00058
00059
00060
00061 PRINT1 (("\nin cholmod_rowfac, kstart %d kend %d stype %d\n",
00062 kstart, kend, A->stype)) ;
00063 DEBUG (CHOLMOD(dump_factor) (L, "Initial L", Common)) ;
00064
00065 n = A->nrow ;
00066 stype = A->stype ;
00067
00068 if (stype > 0)
00069 {
00070
00071 Fp = NULL ;
00072 Fi = NULL ;
00073 Fx = NULL ;
00074 Fz = NULL ;
00075 Fnz = NULL ;
00076 Fpacked = TRUE ;
00077 }
00078 else
00079 {
00080
00081 Fp = F->p ;
00082 Fi = F->i ;
00083 Fx = F->x ;
00084 Fz = F->z ;
00085 Fnz = F->nz ;
00086 Fpacked = F->packed ;
00087 }
00088
00089 Ap = A->p ;
00090 Ai = A->i ;
00091 Ax = A->x ;
00092 Az = A->z ;
00093 Anz = A->nz ;
00094 packed = A->packed ;
00095 sorted = A->sorted ;
00096
00097 use_dbound = IS_GT_ZERO (Common->dbound) ;
00098
00099
00100 is_ll = L->is_ll ;
00101 if (L->xtype == CHOLMOD_PATTERN)
00102 {
00103
00104
00105
00106
00107
00108 CHOLMOD(change_factor) (A->xtype, is_ll, FALSE, FALSE, TRUE, L, Common);
00109 if (Common->status < CHOLMOD_OK)
00110 {
00111
00112 return (FALSE) ;
00113 }
00114 ASSERT (L->minor == (size_t) n) ;
00115 }
00116 else if (kstart == 0 && kend == (size_t) n)
00117 {
00118
00119
00120
00121
00122 L->minor = n ;
00123 Lnz = L->nz ;
00124 for (k = 0 ; k < n ; k++)
00125 {
00126 Lnz [k] = 1 ;
00127 }
00128 }
00129
00130 ASSERT (is_ll == L->is_ll) ;
00131 ASSERT (L->xtype != CHOLMOD_PATTERN) ;
00132 DEBUG (CHOLMOD(dump_factor) (L, "L ready", Common)) ;
00133 DEBUG (CHOLMOD(dump_sparse) (A, "A ready", Common)) ;
00134 DEBUG (if (stype == 0) CHOLMOD(dump_sparse) (F, "F ready", Common)) ;
00135
00136
00137 Lp = L->p ;
00138 ASSERT (Lp != NULL) ;
00139
00140
00141 Lnz = L->nz ;
00142 Lnext = L->next ;
00143 Li = L->i ;
00144 Lx = L->x ;
00145 Lz = L->z ;
00146 nzmax = L->nzmax ;
00147 ASSERT (Lnz != NULL && Li != NULL && Lx != NULL) ;
00148
00149
00150
00151
00152
00153 Iwork = Common->Iwork ;
00154 Stack = Iwork ;
00155 Flag = Common->Flag ;
00156 Wx = Common->Xwork ;
00157
00158 Wz = Wx + n ;
00159 mark = Common->mark ;
00160 ASSERT ((Int) Common->xworksize >= (L->xtype == CHOLMOD_REAL ? 1:2)*n) ;
00161
00162
00163
00164
00165
00166 #ifdef MASK
00167 #define NEXT(k) k = RLinkUp [k]
00168 #else
00169 #define NEXT(k) k++
00170 #endif
00171
00172 for (k = kstart ; k < ((Int) kend) ; NEXT(k))
00173 {
00174 PRINT1 (("\n===============K "ID" Lnz [k] "ID"\n", k, Lnz [k])) ;
00175
00176
00177
00178
00179
00180
00181 ASSERT (Lnz [k] == 1) ;
00182 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 2*n, Common)) ;
00183
00184 top = n ;
00185 Flag [k] = mark ;
00186
00187
00188 #define PARENT(i) (Lnz [i] > 1) ? (Li [Lp [i] + 1]) : EMPTY
00189
00190 if (stype > 0)
00191 {
00192
00193 p = Ap [k] ;
00194 pend = (packed) ? (Ap [k+1]) : (p + Anz [k]) ;
00195
00196 #define SCATTER ASSIGN(Wx,Wz,i, Ax,Az,p)
00197 SUBTREE ;
00198 #undef SCATTER
00199 }
00200 else
00201 {
00202
00203 pf = Fp [k] ;
00204 pfend = (Fpacked) ? (Fp [k+1]) : (pf + Fnz [k]) ;
00205 for ( ; pf < pfend ; pf++)
00206 {
00207
00208 t = Fi [pf] ;
00209
00210 ASSIGN (fx, fz, 0, Fx, Fz, pf) ;
00211 p = Ap [t] ;
00212 pend = (packed) ? (Ap [t+1]) : (p + Anz [t]) ;
00213 multadds = 0 ;
00214
00215 #define SCATTER MULTADD (Wx,Wz,i, Ax,Az,p, fx,fz,0) ; multadds++ ;
00216 SUBTREE ;
00217 #undef SCATTER
00218 #ifdef REAL
00219 fl += 2 * ((double) multadds) ;
00220 #else
00221 fl += 8 * ((double) multadds) ;
00222 #endif
00223 }
00224 }
00225
00226 #undef PARENT
00227
00228
00229
00230
00231
00232 #ifdef MASK
00233
00234 if (mask != NULL)
00235 {
00236
00237 #if 0
00238
00239 for (p = n; p > top;)
00240 {
00241 i = Stack [--p] ;
00242 if ( mask [i] >= 0 )
00243 {
00244 CLEAR (Wx,Wz,i) ;
00245 }
00246 }
00247 #endif
00248
00249 for (s = top ; s < n ; s++)
00250 {
00251 i = Stack [s] ;
00252 if (mask [i] >= 0)
00253 {
00254 CLEAR (Wx,Wz,i) ;
00255 }
00256 }
00257
00258 }
00259 #endif
00260
00261
00262
00263
00264 mark = CHOLMOD(clear_flag) (Common) ;
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283 ADD_REAL (dk,0, Wx,k, beta,0) ;
00284
00285 #ifndef REAL
00286
00287
00288
00289
00290
00291 dk_imaginary = (stype > 0) ? (IMAG_IS_NONZERO (Wx,Wz,k)) : FALSE ;
00292 #endif
00293
00294
00295 CLEAR (Wx,Wz,k) ;
00296
00297 for (s = top ; s < n ; s++)
00298 {
00299
00300 i = Stack [s] ;
00301
00302
00303 ASSIGN (yx,yz,0, Wx,Wz,i) ;
00304
00305
00306 CLEAR (Wx,Wz,i) ;
00307
00308 lnz = Lnz [i] ;
00309 p = Lp [i] ;
00310 ASSERT (lnz > 0 && Li [p] == i) ;
00311 pend = p + lnz ;
00312
00313
00314 ASSIGN_REAL (di,0, Lx,p) ;
00315
00316 if (i >= (Int) L->minor || IS_ZERO (di [0]))
00317 {
00318
00319
00320 CLEAR (lx,lz,0) ;
00321 p = pend ;
00322 }
00323 else if (is_ll)
00324 {
00325 #ifdef REAL
00326 fl += 2 * ((double) (pend - p - 1)) + 3 ;
00327 #else
00328 fl += 8 * ((double) (pend - p - 1)) + 6 ;
00329 #endif
00330
00331
00332
00333 DIV_REAL (yx,yz,0, yx,yz,0, di,0) ;
00334 for (p++ ; p < pend ; p++)
00335 {
00336
00337 MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
00338 }
00339
00340
00341 ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
00342
00343 LLDOT (dk,0, yx,yz,0) ;
00344 }
00345 else
00346 {
00347 #ifdef REAL
00348 fl += 2 * ((double) (pend - p - 1)) + 3 ;
00349 #else
00350 fl += 8 * ((double) (pend - p - 1)) + 6 ;
00351 #endif
00352
00353 for (p++ ; p < pend ; p++)
00354 {
00355
00356 MULTSUB (Wx,Wz,Li[p], Lx,Lz,p, yx,yz,0) ;
00357 }
00358
00359 #ifdef REAL
00360
00361 lx [0] = yx [0] / di [0] ;
00362
00363 dk [0] -= lx [0] * yx [0] ;
00364 #else
00365
00366 ASSIGN_CONJ (lx,lz,0, yx,yz,0) ;
00367
00368 DIV_REAL (lx,lz,0, lx,lz,0, di,0) ;
00369
00370 LDLDOT (dk,0, yx,yz,0, di,0) ;
00371 #endif
00372 }
00373
00374
00375 if (p >= Lp [Lnext [i]])
00376 {
00377
00378 PRINT1 (("Factor Colrealloc "ID", old Lnz "ID"\n", i, Lnz [i]));
00379 if (!CHOLMOD(reallocate_column) (i, lnz + 1, L, Common))
00380 {
00381
00382 for (i = 0 ; i < n ; i++)
00383 {
00384
00385 CLEAR (Wx,Wz,i) ;
00386 }
00387 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
00388 return (FALSE) ;
00389 }
00390 Li = L->i ;
00391 Lx = L->x ;
00392 Lz = L->z ;
00393 p = Lp [i] + lnz ;
00394 ASSERT (p < Lp [Lnext [i]]) ;
00395 }
00396
00397
00398 Li [p] = k ;
00399
00400 ASSIGN (Lx,Lz,p, lx,lz,0) ;
00401 Lnz [i]++ ;
00402 }
00403
00404
00405
00406
00407
00408 p = Lp [k] ;
00409 Li [p] = k ;
00410
00411 if (k >= (Int) L->minor)
00412 {
00413
00414 dk [0] = 0 ;
00415 }
00416 else if (use_dbound)
00417 {
00418
00419 dk [0] = CHOLMOD(dbound) (is_ll ? fabs (dk [0]) : dk [0], Common) ;
00420 }
00421 else if ((is_ll ? (IS_LE_ZERO (dk [0])) : (IS_ZERO (dk [0])))
00422 #ifndef REAL
00423 || dk_imaginary
00424 #endif
00425 )
00426 {
00427
00428 dk [0] = 0 ;
00429 L->minor = k ;
00430 ERROR (CHOLMOD_NOT_POSDEF, "not positive definite") ;
00431 }
00432
00433 if (is_ll)
00434 {
00435
00436 dk [0] = sqrt (dk [0]) ;
00437 }
00438
00439
00440 ASSIGN_REAL (Lx,p, dk,0) ;
00441 CLEAR_IMAG (Lx,Lz,p) ;
00442 }
00443
00444 #undef NEXT
00445
00446 if (is_ll) fl += MAX ((Int) kend - (Int) kstart, 0) ;
00447 Common->rowfacfl = fl ;
00448
00449 DEBUG (CHOLMOD(dump_factor) (L, "final cholmod_rowfac", Common)) ;
00450 ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, n, Common)) ;
00451 return (TRUE) ;
00452 }
00453 #undef PATTERN
00454 #undef REAL
00455 #undef COMPLEX
00456 #undef ZOMPLEX