| Austin Schuh | 9049e20 | 2022-02-20 17:34:16 -0800 | [diff] [blame^] | 1 | #ifndef PARDISO_H |
| 2 | #define PARDISO_H |
| 3 | |
| 4 | #ifdef __cplusplus |
| 5 | extern "C" { |
| 6 | #endif |
| 7 | |
| 8 | #include "lin_alg.h" |
| 9 | #include "kkt.h" |
| 10 | |
| 11 | /** |
| 12 | * Pardiso solver structure |
| 13 | * |
| 14 | * NB: If we use Pardiso, we suppose that EMBEDDED is not enabled |
| 15 | */ |
| 16 | typedef struct pardiso pardiso_solver; |
| 17 | |
| 18 | struct pardiso { |
| 19 | enum linsys_solver_type type; |
| 20 | |
| 21 | /** |
| 22 | * @name Functions |
| 23 | * @{ |
| 24 | */ |
| 25 | c_int (*solve)(struct pardiso * self, c_float * b); |
| 26 | |
| 27 | void (*free)(struct pardiso * self); ///< Free workspace (only if desktop) |
| 28 | |
| 29 | c_int (*update_matrices)(struct pardiso * self, const csc *P, const csc *A); ///< Update solver matrices |
| 30 | c_int (*update_rho_vec)(struct pardiso * self, const c_float * rho_vec); ///< Update rho_vec parameter |
| 31 | |
| 32 | c_int nthreads; |
| 33 | /** @} */ |
| 34 | |
| 35 | |
| 36 | /** |
| 37 | * @name Attributes |
| 38 | * @{ |
| 39 | */ |
| 40 | // Attributes |
| 41 | csc *KKT; ///< KKT matrix (in CSR format!) |
| 42 | c_int *KKT_i; ///< KKT column indices in 1-indexing for Pardiso |
| 43 | c_int *KKT_p; ///< KKT row pointers in 1-indexing for Pardiso |
| 44 | c_float *bp; ///< workspace memory for solves (rhs) |
| 45 | c_float *sol; ///< solution to the KKT system |
| 46 | c_float *rho_inv_vec; ///< parameter vector |
| 47 | c_float sigma; ///< scalar parameter |
| 48 | c_int polish; ///< polishing flag |
| 49 | c_int n; ///< number of QP variables |
| 50 | c_int m; ///< number of QP constraints |
| 51 | |
| 52 | // Pardiso variables |
| 53 | void *pt[64]; ///< internal solver memory pointer pt |
| 54 | c_int iparm[64]; ///< Pardiso control parameters |
| 55 | c_int nKKT; ///< dimension of the linear system |
| 56 | c_int mtype; ///< matrix type (-2 for real and symmetric indefinite) |
| 57 | c_int nrhs; ///< number of right-hand sides (1 for our needs) |
| 58 | c_int maxfct; ///< maximum number of factors (1 for our needs) |
| 59 | c_int mnum; ///< indicates matrix for the solution phase (1 for our needs) |
| 60 | c_int phase; ///< control the execution phases of the solver |
| 61 | c_int error; ///< the error indicator (0 for no error) |
| 62 | c_int msglvl; ///< Message level information (0 for no output) |
| 63 | c_int idum; ///< dummy integer |
| 64 | c_float fdum; ///< dummy float |
| 65 | |
| 66 | // These are required for matrix updates |
| 67 | c_int * Pdiag_idx, Pdiag_n; ///< index and number of diagonal elements in P |
| 68 | c_int * PtoKKT, * AtoKKT; ///< Index of elements from P and A to KKT matrix |
| 69 | c_int * rhotoKKT; ///< Index of rho places in KKT matrix |
| 70 | |
| 71 | /** @} */ |
| 72 | }; |
| 73 | |
| 74 | |
| 75 | /** |
| 76 | * Initialize Pardiso Solver |
| 77 | * |
| 78 | * @param s Pointer to a private structure |
| 79 | * @param P Cost function matrix (upper triangular form) |
| 80 | * @param A Constraints matrix |
| 81 | * @param sigma Algorithm parameter. If polish, then sigma = delta. |
| 82 | * @param rho_vec Algorithm parameter. If polish, then rho_vec = OSQP_NULL. |
| 83 | * @param polish Flag whether we are initializing for polish or not |
| 84 | * @return Exitflag for error (0 if no errors) |
| 85 | */ |
| 86 | 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); |
| 87 | |
| 88 | |
| 89 | /** |
| 90 | * Solve linear system and store result in b |
| 91 | * @param s Linear system solver structure |
| 92 | * @param b Right-hand side |
| 93 | * @return Exitflag |
| 94 | */ |
| 95 | c_int solve_linsys_pardiso(pardiso_solver * s, c_float * b); |
| 96 | |
| 97 | |
| 98 | /** |
| 99 | * Update linear system solver matrices |
| 100 | * @param s Linear system solver structure |
| 101 | * @param P Matrix P |
| 102 | * @param A Matrix A |
| 103 | * @return Exitflag |
| 104 | */ |
| 105 | c_int update_linsys_solver_matrices_pardiso(pardiso_solver * s, const csc *P, const csc *A); |
| 106 | |
| 107 | |
| 108 | /** |
| 109 | * Update rho parameter in linear system solver structure |
| 110 | * @param s Linear system solver structure |
| 111 | * @param rho_vec new rho_vec value |
| 112 | * @return exitflag |
| 113 | */ |
| 114 | c_int update_linsys_solver_rho_vec_pardiso(pardiso_solver * s, const c_float * rho_vec); |
| 115 | |
| 116 | |
| 117 | /** |
| 118 | * Free linear system solver |
| 119 | * @param s linear system solver object |
| 120 | */ |
| 121 | void free_linsys_solver_pardiso(pardiso_solver * s); |
| 122 | |
| 123 | #ifdef __cplusplus |
| 124 | } |
| 125 | #endif |
| 126 | |
| 127 | #endif |