1 /*
2 * basiclu_get_factors.c
3 *
4 * Copyright (C) 2016-2018 ERGO-Code
5 *
6 */
7
8 #include "lu_internal.h"
9
basiclu_get_factors(lu_int istore[],double xstore[],lu_int Li[],double Lx[],lu_int Ui[],double Ux[],lu_int Wi[],double Wx[],lu_int rowperm[],lu_int colperm[],lu_int Lcolptr[],lu_int Lrowidx[],double Lvalue_[],lu_int Ucolptr[],lu_int Urowidx[],double Uvalue_[])10 lu_int basiclu_get_factors
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 lu_int rowperm[],
21 lu_int colperm[],
22 lu_int Lcolptr[],
23 lu_int Lrowidx[],
24 double Lvalue_[],
25 lu_int Ucolptr[],
26 lu_int Urowidx[],
27 double Uvalue_[]
28 )
29 {
30 struct lu this;
31 lu_int m, status;
32
33 status = lu_load(&this, istore, xstore, Li, Lx, Ui, Ux, Wi, Wx);
34 if (status != BASICLU_OK)
35 return status;
36 if (this.nupdate != 0)
37 {
38 status = BASICLU_ERROR_invalid_call;
39 return lu_save(&this, istore, xstore, status);
40 }
41 m = this.m;
42
43 if (rowperm)
44 memcpy(rowperm, this.pivotrow, m*sizeof(lu_int));
45 if (colperm)
46 memcpy(colperm, this.pivotcol, m*sizeof(lu_int));
47
48 if (Lcolptr && Lrowidx && Lvalue_)
49 {
50 const lu_int *Lbegin_p = this.Lbegin_p;
51 const lu_int *Ltbegin_p = this.Ltbegin_p;
52 const lu_int *Lindex = this.Lindex;
53 const double *Lvalue = this.Lvalue;
54 const lu_int *p = this.p;
55 lu_int *colptr = this.iwork1; /* size m workspace */
56 lu_int i, k, put, pos;
57
58 /*
59 * L[:,k] will hold the elimination factors from the k-th pivot step.
60 * First set the column pointers and store the unit diagonal elements
61 * at the front of each column. Then scatter each row of L' into the
62 * columnwise L so that the row indices become sorted.
63 */
64 put = 0;
65 for (k = 0; k < m; k++)
66 {
67 Lcolptr[k] = put;
68 Lrowidx[put] = k;
69 Lvalue_[put++] = 1.0;
70 colptr[p[k]] = put; /* next free position in column */
71 put += Lbegin_p[k+1] - Lbegin_p[k] - 1;
72 /* subtract 1 because internal storage uses (-1) terminators */
73 }
74 Lcolptr[m] = put;
75 assert(put == this.Lnz+m);
76
77 for (k = 0; k < m; k++)
78 {
79 for (pos = Ltbegin_p[k]; (i = Lindex[pos]) >= 0; pos++)
80 {
81 put = colptr[i]++;
82 Lrowidx[put] = k;
83 Lvalue_[put] = Lvalue[pos];
84 }
85 }
86
87 #ifndef NDEBUG
88 for (k = 0; k < m; k++)
89 {
90 assert(colptr[p[k]] == Lcolptr[k+1]);
91 }
92 #endif
93 }
94
95 if (Ucolptr && Urowidx && Uvalue_)
96 {
97 const lu_int *Wbegin = this.Wbegin;
98 const lu_int *Wend = this.Wend;
99 const lu_int *Windex = this.Windex;
100 const double *Wvalue = this.Wvalue;
101 const double *col_pivot = this.col_pivot;
102 const lu_int *pivotcol = this.pivotcol;
103 lu_int *colptr = this.iwork1; /* size m workspace */
104 lu_int j, k, put, pos;
105
106 /*
107 * U[:,k] will hold the column of B from the k-th pivot step.
108 * First set the column pointers and store the pivot element at the end
109 * of each column. Then scatter each row of U' into the columnwise U so
110 * that the row indices become sorted.
111 */
112 memset(colptr, 0, m*sizeof(lu_int)); /* column counts */
113 for (j = 0; j < m; j++)
114 {
115 for (pos = Wbegin[j]; pos < Wend[j]; pos++)
116 colptr[Windex[pos]]++;
117 }
118 put = 0;
119 for (k = 0; k < m; k++) /* set column pointers */
120 {
121 j = pivotcol[k];
122 Ucolptr[k] = put;
123 put += colptr[j];
124 colptr[j] = Ucolptr[k]; /* next free position in column */
125 Urowidx[put] = k;
126 Uvalue_[put++] = col_pivot[j];
127 }
128 Ucolptr[m] = put;
129 assert(put == this.Unz+m);
130 for (k = 0; k < m; k++) /* scatter row k */
131 {
132 j = pivotcol[k];
133 for (pos = Wbegin[j]; pos < Wend[j]; pos++)
134 {
135 put = colptr[Windex[pos]]++;
136 Urowidx[put] = k;
137 Uvalue_[put] = Wvalue[pos];
138 }
139 }
140
141 #ifndef NDEBUG
142 for (k = 0; k < m; k++) assert(colptr[pivotcol[k]] == Ucolptr[k+1]-1);
143 for (k = 0; k < m; k++) assert(Urowidx[Ucolptr[k+1]-1] == k);
144 #endif
145 }
146
147 return BASICLU_OK;
148 }
149