1 #ifndef HEADER_LUSOL 2 #define HEADER_LUSOL 3 4 /* Include necessary libraries */ 5 /* ------------------------------------------------------------------------- */ 6 #include <stdio.h> 7 #include "commonlib.h" 8 9 /* Version information */ 10 /* ------------------------------------------------------------------------- */ 11 #define LUSOL_VERMAJOR 2 12 #define LUSOL_VERMINOR 2 13 #define LUSOL_RELEASE 2 14 #define LUSOL_BUILD 0 15 16 /* Dynamic memory management macros */ 17 /* ------------------------------------------------------------------------- */ 18 #ifdef MATLAB 19 #define LUSOL_MALLOC(bytesize) mxMalloc(bytesize) 20 #define LUSOL_CALLOC(count, recsize) mxCalloc(count, recsize) 21 #define LUSOL_REALLOC(ptr, bytesize) mxRealloc((void *) ptr, bytesize) 22 #define LUSOL_FREE(ptr) {mxFree(ptr); ptr=NULL;} 23 #else 24 #define LUSOL_MALLOC(bytesize) malloc(bytesize) 25 #define LUSOL_CALLOC(count, recsize) calloc(count, recsize) 26 #define LUSOL_REALLOC(ptr, bytesize) realloc((void *) ptr, bytesize) 27 #define LUSOL_FREE(ptr) {free(ptr); ptr=NULL;} 28 #endif 29 30 /* Performance compiler options */ 31 /* ------------------------------------------------------------------------- */ 32 #if 1 33 #define ForceInitialization /* Satisfy compilers, check during debugging! */ 34 #define LUSOLFastDenseIndex /* Increment the linearized dense address */ 35 #define LUSOLFastClear /* Use intrinsic functions for memory zeroing */ 36 #define LUSOLFastMove /* Use intrinsic functions for memory moves */ 37 #define LUSOLFastCopy /* Use intrinsic functions for memory copy */ 38 #define LUSOLFastSolve /* Use pointer operations in equation solving */ 39 #define LUSOLSafeFastUpdate /* Use separate array for LU6L result storage */ 40 /*#define UseOld_LU6CHK_20040510 */ 41 /*#define AlwaysSeparateHamaxR */ /* Enabled when the pivot model is fixed */ 42 #if 0 43 #define ForceRowBasedL0 /* Create a row-sorted version of L0 */ 44 #endif 45 /* #define SetSmallToZero*/ 46 /* #define DoTraceL0 */ 47 #endif 48 /*#define UseTimer */ 49 50 51 /* Legacy compatibility and testing options (Fortran-LUSOL) */ 52 /* ------------------------------------------------------------------------- */ 53 #if 0 54 #define LegacyTesting 55 #define StaticMemAlloc /* Preallocated vs. dynamic memory allocation */ 56 #define ClassicdiagU /* Store diagU at end of a */ 57 #define ClassicHamaxR /* Store H+AmaxR at end of a/indc/indr */ 58 #endif 59 60 61 /* General constants and data type definitions */ 62 /* ------------------------------------------------------------------------- */ 63 #define LUSOL_ARRAYOFFSET 1 64 #ifndef ZERO 65 #define ZERO 0 66 #endif 67 #ifndef ONE 68 #define ONE 1 69 #endif 70 #ifndef FALSE 71 #define FALSE 0 72 #endif 73 #ifndef TRUE 74 #define TRUE 1 75 #endif 76 #ifndef NULL 77 #define NULL 0 78 #endif 79 #ifndef REAL 80 #define REAL double 81 #endif 82 #ifndef REALXP 83 #define REALXP long double 84 #endif 85 #ifndef MYBOOL 86 #define MYBOOL unsigned char 87 #endif 88 89 90 /* User-settable default parameter values */ 91 /* ------------------------------------------------------------------------- */ 92 #define LUSOL_DEFAULT_GAMMA 2.0 93 #define LUSOL_SMALLNUM 1.0e-20 /* IAEE doubles have precision 2.22e-16 */ 94 #define LUSOL_BIGNUM 1.0e+20 95 #define LUSOL_MINDELTA_FACTOR 4 96 #define LUSOL_MINDELTA_a 10000 97 #if 1 98 #define LUSOL_MULT_nz_a 2 /* Suggested by Yin Zhang */ 99 #else 100 #define LUSOL_MULT_nz_a 5 /* Could consider 6 or 7 */ 101 #endif 102 #define LUSOL_MINDELTA_rc 1000 103 #define LUSOL_DEFAULT_SMARTRATIO 0.667 104 105 /* Fixed system parameters (changeable only by developers) */ 106 /* ------------------------------------------------------------------------- */ 107 108 /* parmlu INPUT parameters: */ 109 #define LUSOL_RP_SMARTRATIO 0 110 #define LUSOL_RP_FACTORMAX_Lij 1 111 #define LUSOL_RP_UPDATEMAX_Lij 2 112 #define LUSOL_RP_ZEROTOLERANCE 3 113 #define LUSOL_RP_SMALLDIAG_U 4 114 #define LUSOL_RP_EPSDIAG_U 5 115 #define LUSOL_RP_COMPSPACE_U 6 116 #define LUSOL_RP_MARKOWITZ_CONLY 7 117 #define LUSOL_RP_MARKOWITZ_DENSE 8 118 #define LUSOL_RP_GAMMA 9 119 120 /* parmlu OUPUT parameters: */ 121 #define LUSOL_RP_MAXELEM_A 10 122 #define LUSOL_RP_MAXMULT_L 11 123 #define LUSOL_RP_MAXELEM_U 12 124 #define LUSOL_RP_MAXELEM_DIAGU 13 125 #define LUSOL_RP_MINELEM_DIAGU 14 126 #define LUSOL_RP_MAXELEM_TCP 15 127 #define LUSOL_RP_GROWTHRATE 16 128 #define LUSOL_RP_USERDATA_1 17 129 #define LUSOL_RP_USERDATA_2 18 130 #define LUSOL_RP_USERDATA_3 19 131 #define LUSOL_RP_RESIDUAL_U 20 132 #define LUSOL_RP_LASTITEM LUSOL_RP_RESIDUAL_U 133 134 /* luparm INPUT parameters: */ 135 #define LUSOL_IP_USERDATA_0 0 136 #define LUSOL_IP_PRINTUNIT 1 137 #define LUSOL_IP_PRINTLEVEL 2 138 #define LUSOL_IP_MARKOWITZ_MAXCOL 3 139 #define LUSOL_IP_SCALAR_NZA 4 140 #define LUSOL_IP_UPDATELIMIT 5 141 #define LUSOL_IP_PIVOTTYPE 6 142 #define LUSOL_IP_ACCELERATION 7 143 #define LUSOL_IP_KEEPLU 8 144 #define LUSOL_IP_SINGULARLISTSIZE 9 145 146 /* luparm OUTPUT parameters: */ 147 #define LUSOL_IP_INFORM 10 148 #define LUSOL_IP_SINGULARITIES 11 149 #define LUSOL_IP_SINGULARINDEX 12 150 #define LUSOL_IP_MINIMUMLENA 13 151 #define LUSOL_IP_MAXLEN 14 152 #define LUSOL_IP_UPDATECOUNT 15 153 #define LUSOL_IP_RANK_U 16 154 #define LUSOL_IP_COLCOUNT_DENSE1 17 155 #define LUSOL_IP_COLCOUNT_DENSE2 18 156 #define LUSOL_IP_COLINDEX_DUMIN 19 157 #define LUSOL_IP_COLCOUNT_L0 20 158 #define LUSOL_IP_NONZEROS_L0 21 159 #define LUSOL_IP_NONZEROS_U0 22 160 #define LUSOL_IP_NONZEROS_L 23 161 #define LUSOL_IP_NONZEROS_U 24 162 #define LUSOL_IP_NONZEROS_ROW 25 163 #define LUSOL_IP_COMPRESSIONS_LU 26 164 #define LUSOL_IP_MARKOWITZ_MERIT 27 165 #define LUSOL_IP_TRIANGROWS_U 28 166 #define LUSOL_IP_TRIANGROWS_L 29 167 #define LUSOL_IP_FTRANCOUNT 30 168 #define LUSOL_IP_BTRANCOUNT 31 169 #define LUSOL_IP_ROWCOUNT_L0 32 170 #define LUSOL_IP_LASTITEM LUSOL_IP_ROWCOUNT_L0 171 172 173 /* Macros for matrix-based access for dense part of A and timer mapping */ 174 /* ------------------------------------------------------------------------- */ 175 #define DAPOS(row, col) (row + (col-1)*LDA) 176 #define timer(text, id) LUSOL_timer(LUSOL, id, text) 177 178 179 /* Parameter/option defines */ 180 /* ------------------------------------------------------------------------- */ 181 #define LUSOL_MSG_NONE -1 182 #define LUSOL_MSG_SINGULARITY 0 183 #define LUSOL_MSG_STATISTICS 10 184 #define LUSOL_MSG_PIVOT 50 185 186 #define LUSOL_BASEORDER 0 187 #define LUSOL_OTHERORDER 1 188 #define LUSOL_AUTOORDER 2 189 #define LUSOL_ACCELERATE_L0 4 190 #define LUSOL_ACCELERATE_U 8 191 192 #define LUSOL_PIVMOD_NOCHANGE -2 /* Don't change active pivoting model */ 193 #define LUSOL_PIVMOD_DEFAULT -1 /* Set pivoting model to default */ 194 #define LUSOL_PIVMOD_TPP 0 /* Threshold Partial pivoting (normal) */ 195 #define LUSOL_PIVMOD_TRP 1 /* Threshold Rook pivoting */ 196 #define LUSOL_PIVMOD_TCP 2 /* Threshold Complete pivoting */ 197 #define LUSOL_PIVMOD_TSP 3 /* Threshold Symmetric pivoting */ 198 #define LUSOL_PIVMOD_MAX LUSOL_PIVMOD_TSP 199 200 #define LUSOL_PIVTOL_NOCHANGE 0 201 #define LUSOL_PIVTOL_BAGGY 1 202 #define LUSOL_PIVTOL_LOOSE 2 203 #define LUSOL_PIVTOL_NORMAL 3 204 #define LUSOL_PIVTOL_SLIM 4 205 #define LUSOL_PIVTOL_TIGHT 5 206 #define LUSOL_PIVTOL_SUPER 6 207 #define LUSOL_PIVTOL_CORSET 7 208 #define LUSOL_PIVTOL_DEFAULT LUSOL_PIVTOL_SLIM 209 #define LUSOL_PIVTOL_MAX LUSOL_PIVTOL_CORSET 210 211 #define LUSOL_UPDATE_OLDEMPTY 0 /* No/empty current column. */ 212 #define LUSOL_UPDATE_OLDNONEMPTY 1 /* Current column need not have been empty. */ 213 #define LUSOL_UPDATE_NEWEMPTY 0 /* New column is taken to be zero. */ 214 #define LUSOL_UPDATE_NEWNONEMPTY 1 /* v(*) contains the new column; 215 on exit, v(*) satisfies L*v = a(new). */ 216 #define LUSOL_UPDATE_USEPREPARED 2 /* v(*) must satisfy L*v = a(new). */ 217 218 #define LUSOL_SOLVE_Lv_v 1 /* v solves L v = v(input). w is not touched. */ 219 #define LUSOL_SOLVE_Ltv_v 2 /* v solves L'v = v(input). w is not touched. */ 220 #define LUSOL_SOLVE_Uw_v 3 /* w solves U w = v. v is not altered. */ 221 #define LUSOL_SOLVE_Utv_w 4 /* v solves U'v = w. w is destroyed. */ 222 #define LUSOL_SOLVE_Aw_v 5 /* w solves A w = v. v is altered as in 1. */ 223 #define LUSOL_FTRAN LUSOL_SOLVE_Aw_v 224 #define LUSOL_SOLVE_Atv_w 6 /* v solves A'v = w. w is destroyed. */ 225 #define LUSOL_BTRAN LUSOL_SOLVE_Atv_w 226 227 /* If mode = 3,4,5,6, v and w must not be the same arrays. 228 If lu1fac has just been used to factorize a symmetric matrix A 229 (which must be definite or quasi-definite), the factors A = L U 230 may be regarded as A = LDL', where D = diag(U). In such cases, 231 the following (faster) solve codes may be used: */ 232 #define LUSOL_SOLVE_Av_v 7 /* v solves A v = L D L'v = v(input). w is not touched. */ 233 #define LUSOL_SOLVE_LDLtv_v 8 /* v solves L |D| L'v = v(input). w is not touched. */ 234 235 #define LUSOL_INFORM_RANKLOSS -1 236 #define LUSOL_INFORM_LUSUCCESS 0 237 #define LUSOL_INFORM_LUSINGULAR 1 238 #define LUSOL_INFORM_LUUNSTABLE 2 239 #define LUSOL_INFORM_ADIMERR 3 240 #define LUSOL_INFORM_ADUPLICATE 4 241 #define LUSOL_INFORM_ANEEDMEM 7 /* Set lena >= luparm[LUSOL_IP_MINIMUMLENA] */ 242 #define LUSOL_INFORM_FATALERR 8 243 #define LUSOL_INFORM_NOPIVOT 9 /* No diagonal pivot found with TSP or TDP. */ 244 #define LUSOL_INFORM_NOMEMLEFT 10 245 246 #define LUSOL_INFORM_MIN LUSOL_INFORM_RANKLOSS 247 #define LUSOL_INFORM_MAX LUSOL_INFORM_NOMEMLEFT 248 249 #define LUSOL_INFORM_GETLAST 10 /* Code for LUSOL_informstr. */ 250 #define LUSOL_INFORM_SERIOUS LUSOL_INFORM_LUUNSTABLE 251 252 253 /* Prototypes for call-back functions */ 254 /* ------------------------------------------------------------------------- */ 255 typedef void LUSOLlogfunc(void *lp, void *userhandle, char *buf); 256 257 258 /* Sparse matrix data */ 259 typedef struct _LUSOLmat { 260 REAL *a; 261 int *lenx, *indr, *indc, *indx; 262 } LUSOLmat; 263 264 265 /* The main LUSOL data record */ 266 /* ------------------------------------------------------------------------- */ 267 typedef struct _LUSOLrec { 268 269 /* General data */ 270 FILE *outstream; /* Output stream, initialized to STDOUT */ 271 LUSOLlogfunc *writelog; 272 void *loghandle; 273 LUSOLlogfunc *debuginfo; 274 275 /* Parameter storage arrays */ 276 int luparm[LUSOL_IP_LASTITEM + 1]; 277 REAL parmlu[LUSOL_RP_LASTITEM + 1]; 278 279 /* Arrays of length lena+1 */ 280 int lena, nelem; 281 int *indc, *indr; 282 REAL *a; 283 284 /* Arrays of length maxm+1 (row storage) */ 285 int maxm, m; 286 int *lenr, *ip, *iqloc, *ipinv, *locr; 287 288 /* Arrays of length maxn+1 (column storage) */ 289 int maxn, n; 290 int *lenc, *iq, *iploc, *iqinv, *locc; 291 REAL *w, *vLU6L; 292 293 /* List of singular columns, with dynamic size allocation */ 294 int *isingular; 295 296 /* Extra arrays of length n for TCP and keepLU == FALSE */ 297 REAL *Ha, *diagU; 298 int *Hj, *Hk; 299 300 /* Extra arrays of length m for TRP*/ 301 REAL *amaxr; 302 303 /* Extra array for L0 and U stored by row/column for faster btran/ftran */ 304 LUSOLmat *L0; 305 LUSOLmat *U; 306 307 /* Miscellaneous data */ 308 int expanded_a; 309 int replaced_c; 310 int replaced_r; 311 312 } LUSOLrec; 313 314 315 LUSOLrec *LUSOL_create(FILE *outstream, int msgfil, int pivotmodel, int updatelimit); 316 MYBOOL LUSOL_sizeto(LUSOLrec *LUSOL, int init_r, int init_c, int init_a); 317 MYBOOL LUSOL_assign(LUSOLrec *LUSOL, int iA[], int jA[], REAL Aij[], 318 int nzcount, MYBOOL istriplet); 319 void LUSOL_clear(LUSOLrec *LUSOL, MYBOOL nzonly); 320 void LUSOL_free(LUSOLrec *LUSOL); 321 322 LUSOLmat *LUSOL_matcreate(int dim, int nz); 323 void LUSOL_matfree(LUSOLmat **mat); 324 325 int LUSOL_loadColumn(LUSOLrec *LUSOL, int iA[], int jA, REAL Aij[], int nzcount, int offset1); 326 void LUSOL_setpivotmodel(LUSOLrec *LUSOL, int pivotmodel, int initlevel); 327 int LUSOL_factorize(LUSOLrec *LUSOL); 328 int LUSOL_replaceColumn(LUSOLrec *LUSOL, int jcol, REAL v[]); 329 330 MYBOOL LUSOL_tightenpivot(LUSOLrec *LUSOL); 331 MYBOOL LUSOL_addSingularity(LUSOLrec *LUSOL, int singcol, int *inform); 332 int LUSOL_getSingularity(LUSOLrec *LUSOL, int singitem); 333 int LUSOL_findSingularityPosition(LUSOLrec *LUSOL, int singcol); 334 335 char *LUSOL_pivotLabel(LUSOLrec *LUSOL); 336 char *LUSOL_informstr(LUSOLrec *LUSOL, int inform); 337 REAL LUSOL_vecdensity(LUSOLrec *LUSOL, REAL V[]); 338 void LUSOL_report(LUSOLrec *LUSOL, int msglevel, char *format, ...); 339 void LUSOL_timer(LUSOLrec *LUSOL, int timerid, char *text); 340 341 int LUSOL_ftran(LUSOLrec *LUSOL, REAL b[], int NZidx[], MYBOOL prepareupdate); 342 int LUSOL_btran(LUSOLrec *LUSOL, REAL b[], int NZidx[]); 343 344 void LU1FAC(LUSOLrec *LUSOL, int *INFORM); 345 MYBOOL LU1L0(LUSOLrec *LUSOL, LUSOLmat **mat, int *inform); 346 void LU6SOL(LUSOLrec *LUSOL, int MODE, REAL V[], REAL W[], int NZidx[], int *INFORM); 347 void LU8RPC(LUSOLrec *LUSOL, int MODE1, int MODE2, 348 int JREP, REAL V[], REAL W[], 349 int *INFORM, REAL *DIAG, REAL *VNORM); 350 351 void LUSOL_dump(FILE *output, LUSOLrec *LUSOL); 352 353 354 void print_L0(LUSOLrec *LUSOL); 355 356 357 #endif /* HEADER_LUSOL */ 358