1 /*
2 * basiclu_factorize.c
3 *
4 * Copyright (C) 2016-2018 ERGO-Code
5 *
6 */
7
8 #include "lu_internal.h"
9
basiclu_factorize(lu_int istore[],double xstore[],lu_int Li[],double Lx[],lu_int Ui[],double Ux[],lu_int Wi[],double Wx[],const lu_int Bbegin[],const lu_int Bend[],const lu_int Bi[],const double Bx[],lu_int c0ntinue)10 lu_int basiclu_factorize
11 (
12 lu_int istore[],
13 double xstore[],
14 lu_int Li[],
15 double Lx[],
16 lu_int Ui[],
17 double Ux[],
18 lu_int Wi[],
19 double Wx[],
20 const lu_int Bbegin[],
21 const lu_int Bend[],
22 const lu_int Bi[],
23 const double Bx[],
24 lu_int c0ntinue
25 )
26 {
27 struct lu this;
28 lu_int status;
29 double tic[2], elapsed, factor_cost;
30
31 status = lu_load(&this, istore, xstore, Li, Lx, Ui, Ux, Wi, Wx);
32 if (status != BASICLU_OK)
33 return status;
34
35 if (! (Li && Lx && Ui && Ux && Wi && Wx && Bbegin && Bend && Bi && Bx))
36 {
37 status = BASICLU_ERROR_argument_missing;
38 return lu_save(&this, istore, xstore, status);
39 }
40 if (!c0ntinue)
41 {
42 lu_reset(&this);
43 this.task = SINGLETONS;
44 }
45
46 /* continue factorization */
47 switch (this.task)
48 {
49 case SINGLETONS: goto singletons;
50 case SETUP_BUMP: goto setup_bump;
51 case FACTORIZE_BUMP: goto factorize_bump;
52 case BUILD_FACTORS: goto build_factors;
53 }
54 status = BASICLU_ERROR_invalid_call;
55 return lu_save(&this, istore, xstore, status);
56
57 /*
58 * Each of the following four parts of the factorization calls a routine
59 * lu_do_something() which may request reallocation. In this case return
60 * to the caller immediately, keeping the entry point in this.task.
61 */
62
63 singletons:
64 this.task = SINGLETONS;
65 status = lu_singletons(&this, Bbegin, Bend, Bi, Bx);
66 if (status != BASICLU_OK)
67 goto return_to_caller;
68
69 setup_bump:
70 this.task = SETUP_BUMP;
71 status = lu_setup_bump(&this, Bbegin, Bend, Bi, Bx);
72 if (status != BASICLU_OK)
73 goto return_to_caller;
74
75 factorize_bump:
76 this.task = FACTORIZE_BUMP;
77 status = lu_factorize_bump(&this);
78 if (status != BASICLU_OK)
79 goto return_to_caller;
80
81 build_factors:
82 this.task = BUILD_FACTORS;
83 status = lu_build_factors(&this);
84 if (status != BASICLU_OK)
85 goto return_to_caller;
86
87 /* factorization successfully finished */
88 this.task = NO_TASK;
89 this.nupdate = 0; /* make factorization valid */
90 this.ftran_for_update = -1;
91 this.btran_for_update = -1;
92 this.nfactorize++;
93
94 this.condestL = lu_condest(this.m, this.Lbegin, this.Lindex, this.Lvalue,
95 NULL, this.p, 0, this.work1,
96 &this.normL, &this.normestLinv);
97 this.condestU = lu_condest(this.m, this.Ubegin, this.Uindex, this.Uvalue,
98 this.row_pivot, this.p, 1, this.work1,
99 &this.normU, &this.normestUinv);
100
101 /* measure numerical stability of the factorization */
102 lu_residual_test(&this, Bbegin, Bend, Bi, Bx);
103
104 /*
105 * factor_cost is a deterministic measure of the factorization cost.
106 * The parameters have been adjusted such that (on my computer)
107 * 1e-6 * factor_cost =~ time_factorize.
108 *
109 * update_cost measures the accumulated cost of updates/solves compared
110 * to the last factorization. It is computed from
111 *
112 * update_cost = update_cost_numer / update_cost_denom.
113 *
114 * update_cost_denom is fixed here.
115 * update_cost_numer is zero here and increased by solves/updates.
116 *
117 */
118 factor_cost =
119 0.04 * this.m +
120 0.07 * this.matrix_nz +
121 0.20 * this.bump_nz +
122 0.20 * this.nsearch_pivot +
123 0.008 * this.factor_flops;
124
125 this.update_cost_denom = factor_cost * 250;
126
127 if (this.rank < this.m)
128 status = BASICLU_WARNING_singular_matrix;
129
130 return_to_caller:
131 return lu_save(&this, istore, xstore, status);
132 }
133