1 #include "pardiso_interface.h"
2
3 #define MKL_INT c_int
4
5 // Single Dynamic library interface
6 #define MKL_INTERFACE_LP64 0x0
7 #define MKL_INTERFACE_ILP64 0x1
8
9 // Solver Phases
10 #define PARDISO_SYMBOLIC (11)
11 #define PARDISO_NUMERIC (22)
12 #define PARDISO_SOLVE (33)
13 #define PARDISO_CLEANUP (-1)
14
15
16 // Prototypes for Pardiso functions
17 void pardiso(void**, // pt
18 const c_int*, // maxfct
19 const c_int*, // mnum
20 const c_int*, // mtype
21 const c_int*, // phase
22 const c_int*, // n
23 const c_float*, // a
24 const c_int*, // ia
25 const c_int*, // ja
26 c_int*, // perm
27 const c_int*, //nrhs
28 c_int*, // iparam
29 const c_int*, //msglvl
30 c_float*, // b
31 c_float*, // x
32 c_int* // error
33 );
34 c_int mkl_set_interface_layer(c_int);
35 c_int mkl_get_max_threads();
36
37 // Free LDL Factorization structure
free_linsys_solver_pardiso(pardiso_solver * s)38 void free_linsys_solver_pardiso(pardiso_solver *s) {
39 if (s) {
40
41 // Free pardiso solver using internal function
42 s->phase = PARDISO_CLEANUP;
43 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
44 &(s->nKKT), &(s->fdum), s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
45 s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
46
47 if ( s->error != 0 ){
48 #ifdef PRINTING
49 c_eprint("Error during MKL Pardiso cleanup: %d", (int)s->error);
50 #endif
51 }
52 // Check each attribute of the structure and free it if it exists
53 if (s->KKT) csc_spfree(s->KKT);
54 if (s->KKT_i) c_free(s->KKT_i);
55 if (s->KKT_p) c_free(s->KKT_p);
56 if (s->bp) c_free(s->bp);
57 if (s->sol) c_free(s->sol);
58 if (s->rho_inv_vec) c_free(s->rho_inv_vec);
59
60 // These are required for matrix updates
61 if (s->Pdiag_idx) c_free(s->Pdiag_idx);
62 if (s->PtoKKT) c_free(s->PtoKKT);
63 if (s->AtoKKT) c_free(s->AtoKKT);
64 if (s->rhotoKKT) c_free(s->rhotoKKT);
65
66 c_free(s);
67
68 }
69 }
70
71
72 // Initialize factorization structure
init_linsys_solver_pardiso(pardiso_solver ** sp,const csc * P,const csc * A,c_float sigma,const c_float * rho_vec,c_int polish)73 c_int init_linsys_solver_pardiso(pardiso_solver ** sp, const csc * P, const csc * A, c_float sigma, const c_float * rho_vec, c_int polish){
74 c_int i; // loop counter
75 c_int nnzKKT; // Number of nonzeros in KKT
76 // Define Variables
77 c_int n_plus_m; // n_plus_m dimension
78
79
80 // Allocate private structure to store KKT factorization
81 pardiso_solver *s;
82 s = c_calloc(1, sizeof(pardiso_solver));
83 *sp = s;
84
85 // Size of KKT
86 s->n = P->n;
87 s->m = A->m;
88 n_plus_m = s->n + s->m;
89 s->nKKT = n_plus_m;
90
91 // Sigma parameter
92 s->sigma = sigma;
93
94 // Polishing flag
95 s->polish = polish;
96
97 // Link Functions
98 s->solve = &solve_linsys_pardiso;
99 s->free = &free_linsys_solver_pardiso;
100 s->update_matrices = &update_linsys_solver_matrices_pardiso;
101 s->update_rho_vec = &update_linsys_solver_rho_vec_pardiso;
102
103 // Assign type
104 s->type = MKL_PARDISO_SOLVER;
105
106 // Working vector
107 s->bp = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
108
109 // Solution vector
110 s->sol = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
111
112 // Parameter vector
113 s->rho_inv_vec = (c_float *)c_malloc(sizeof(c_float) * n_plus_m);
114
115 // Form KKT matrix
116 if (polish){ // Called from polish()
117 // Use s->rho_inv_vec for storing param2 = vec(delta)
118 for (i = 0; i < A->m; i++){
119 s->rho_inv_vec[i] = sigma;
120 }
121
122 s->KKT = form_KKT(P, A, 1, sigma, s->rho_inv_vec, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL, OSQP_NULL);
123 }
124 else { // Called from ADMM algorithm
125
126 // Allocate vectors of indices
127 s->PtoKKT = c_malloc((P->p[P->n]) * sizeof(c_int));
128 s->AtoKKT = c_malloc((A->p[A->n]) * sizeof(c_int));
129 s->rhotoKKT = c_malloc((A->m) * sizeof(c_int));
130
131 // Use s->rho_inv_vec for storing param2 = rho_inv_vec
132 for (i = 0; i < A->m; i++){
133 s->rho_inv_vec[i] = 1. / rho_vec[i];
134 }
135
136 s->KKT = form_KKT(P, A, 1, sigma, s->rho_inv_vec,
137 s->PtoKKT, s->AtoKKT,
138 &(s->Pdiag_idx), &(s->Pdiag_n), s->rhotoKKT);
139 }
140
141 // Check if matrix has been created
142 if (!(s->KKT)) {
143 #ifdef PRINTING
144 c_eprint("Error in forming KKT matrix");
145 #endif
146 free_linsys_solver_pardiso(s);
147 return OSQP_LINSYS_SOLVER_INIT_ERROR;
148 } else {
149 // Adjust indexing for Pardiso
150 nnzKKT = s->KKT->p[s->KKT->m];
151 s->KKT_i = c_malloc((nnzKKT) * sizeof(c_int));
152 s->KKT_p = c_malloc((s->KKT->m + 1) * sizeof(c_int));
153
154 for(i = 0; i < nnzKKT; i++){
155 s->KKT_i[i] = s->KKT->i[i] + 1;
156 }
157 for(i = 0; i < n_plus_m+1; i++){
158 s->KKT_p[i] = s->KKT->p[i] + 1;
159 }
160
161 }
162
163 // Set MKL interface layer (Long integers if activated)
164 #ifdef DLONG
165 mkl_set_interface_layer(MKL_INTERFACE_ILP64);
166 #else
167 mkl_set_interface_layer(MKL_INTERFACE_LP64);
168 #endif
169
170 // Set Pardiso variables
171 s->mtype = -2; // Real symmetric indefinite matrix
172 s->nrhs = 1; // Number of right hand sides
173 s->maxfct = 1; // Maximum number of numerical factorizations
174 s->mnum = 1; // Which factorization to use
175 s->msglvl = 0; // Do not print statistical information
176 s->error = 0; // Initialize error flag
177 for ( i = 0; i < 64; i++ ) {
178 s->iparm[i] = 0; // Setup Pardiso control parameters
179 s->pt[i] = 0; // Initialize the internal solver memory pointer
180 }
181 s->iparm[0] = 1; // No solver default
182 s->iparm[1] = 3; // Fill-in reordering from OpenMP
183 if (polish) {
184 s->iparm[5] = 1; // Write solution into b
185 } else {
186 s->iparm[5] = 0; // Do NOT write solution into b
187 }
188 /* s->iparm[7] = 2; // Max number of iterative refinement steps */
189 s->iparm[7] = 0; // Number of iterative refinement steps (auto, performs them only if perturbed pivots are obtained)
190 s->iparm[9] = 13; // Perturb the pivot elements with 1E-13
191 s->iparm[34] = 0; // Use Fortran-style indexing for indices
192 /* s->iparm[34] = 1; // Use C-style indexing for indices */
193
194 // Print number of threads
195 s->nthreads = mkl_get_max_threads();
196
197 // Reordering and symbolic factorization
198 s->phase = PARDISO_SYMBOLIC;
199 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
200 &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
201 s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
202 if ( s->error != 0 ){
203 #ifdef PRINTING
204 c_eprint("Error during symbolic factorization: %d", (int)s->error);
205 #endif
206 free_linsys_solver_pardiso(s);
207 *sp = OSQP_NULL;
208 return OSQP_LINSYS_SOLVER_INIT_ERROR;
209 }
210
211 // Numerical factorization
212 s->phase = PARDISO_NUMERIC;
213 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
214 &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
215 s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
216 if ( s->error ){
217 #ifdef PRINTING
218 c_eprint("Error during numerical factorization: %d", (int)s->error);
219 #endif
220 free_linsys_solver_pardiso(s);
221 *sp = OSQP_NULL;
222 return OSQP_LINSYS_SOLVER_INIT_ERROR;
223 }
224
225
226 // No error
227 return 0;
228 }
229
230 // Returns solution to linear system Ax = b with solution stored in b
solve_linsys_pardiso(pardiso_solver * s,c_float * b)231 c_int solve_linsys_pardiso(pardiso_solver * s, c_float * b) {
232 c_int j;
233
234 // Back substitution and iterative refinement
235 s->phase = PARDISO_SOLVE;
236 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
237 &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
238 s->iparm, &(s->msglvl), b, s->sol, &(s->error));
239 if ( s->error != 0 ){
240 #ifdef PRINTING
241 c_eprint("Error during linear system solution: %d", (int)s->error);
242 #endif
243 return 1;
244 }
245
246 if (!(s->polish)) {
247 /* copy x_tilde from s->sol */
248 for (j = 0 ; j < s->n ; j++) {
249 b[j] = s->sol[j];
250 }
251
252 /* compute z_tilde from b and s->sol */
253 for (j = 0 ; j < s->m ; j++) {
254 b[j + s->n] += s->rho_inv_vec[j] * s->sol[j + s->n];
255 }
256 }
257
258 return 0;
259 }
260
261 // Update solver structure with new P and A
update_linsys_solver_matrices_pardiso(pardiso_solver * s,const csc * P,const csc * A)262 c_int update_linsys_solver_matrices_pardiso(pardiso_solver * s, const csc *P, const csc *A) {
263
264 // Update KKT matrix with new P
265 update_KKT_P(s->KKT, P, s->PtoKKT, s->sigma, s->Pdiag_idx, s->Pdiag_n);
266
267 // Update KKT matrix with new A
268 update_KKT_A(s->KKT, A, s->AtoKKT);
269
270 // Perform numerical factorization
271 s->phase = PARDISO_NUMERIC;
272 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
273 &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
274 s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
275
276 // Return exit flag
277 return s->error;
278 }
279
280
update_linsys_solver_rho_vec_pardiso(pardiso_solver * s,const c_float * rho_vec)281 c_int update_linsys_solver_rho_vec_pardiso(pardiso_solver * s, const c_float * rho_vec) {
282 c_int i;
283
284 // Update internal rho_inv_vec
285 for (i = 0; i < s->m; i++){
286 s->rho_inv_vec[i] = 1. / rho_vec[i];
287 }
288
289 // Update KKT matrix with new rho_vec
290 update_KKT_param2(s->KKT, s->rho_inv_vec, s->rhotoKKT, s->m);
291
292 // Perform numerical factorization
293 s->phase = PARDISO_NUMERIC;
294 pardiso (s->pt, &(s->maxfct), &(s->mnum), &(s->mtype), &(s->phase),
295 &(s->nKKT), s->KKT->x, s->KKT_p, s->KKT_i, &(s->idum), &(s->nrhs),
296 s->iparm, &(s->msglvl), &(s->fdum), &(s->fdum), &(s->error));
297
298 // Return exit flag
299 return s->error;
300 }
301