1
2 /* This program solves a sparse linear system Ax = b using LUSOL. */
3
4 #include <stdlib.h>
5 #include <stdio.h>
6 #include <math.h>
7 #include <string.h>
8 #include <time.h>
9 #include "commonlib.h"
10 #include "myblas.h"
11 #include "lusol.h"
12 #include "lusolio.h"
13 #include "lusolmain.h"
14
15 #if (defined WIN32) || (defined WIN64)
_strupr_(char * s)16 void _strupr_(char *s)
17 {
18 _strupr(s);
19 }
20 #else
21 /* Yin Zhang noted that _strupr() is not available on many Unix platforms */
_strupr_(char * s)22 void _strupr_(char *s)
23 {
24 int i;
25 char c;
26 for (i = 0; i < strlen(s); i++) {
27 c = s[i];
28 if (c <= 'z' && c >= 'a') {
29 s[i] = c - 'a' + 'A';
30 }
31 }
32 }
33 #endif
34
getFileName(char * filename,char * test)35 MYBOOL getFileName(char *filename, char *test)
36 {
37 MYBOOL status;
38 status = (MYBOOL) (('-' != test[0]) && (strlen(test) > 1));
39 if(status)
40 strcpy(filename, test);
41 return(status);
42 }
isNum(char val)43 MYBOOL isNum(char val)
44 {
45 int ord;
46 ord = (int) val - 48;
47 return( (MYBOOL) ((ord >= 0) && (ord <= 9)) );
48 }
49
main(int argc,char * argv[],char * envp[])50 void main( int argc, char *argv[], char *envp[] )
51 {
52 /* Output device */
53 FILE *outunit = stdout;
54
55 /* Overall dimensions allocated */
56 int maxm = MAXROWS, maxn = MAXCOLS, maxnz = MAXNZ,
57 replace = 0, randcol = 0;
58 MYBOOL ftran = TRUE;
59
60 /* Storage for A, b */
61 REAL *Aij = NULL, *b = NULL, *xexact = NULL;
62 int *iA = NULL, *jA = NULL;
63
64 /* Storage for LUSOL */
65 LUSOLrec *LUSOL = NULL;
66
67 /* Define local storage variables */
68 int i , inform, j , k , i1 ,
69 m , useropt, lenb, lenx,
70 n , nelem , nnzero;
71 REAL Amax , test ,
72 bnorm , rnorm , xnorm,
73 *rhs = NULL, *r = NULL, *x = NULL;
74 char fileext[50], filename[255], blasname[255];
75 MYBOOL printsolution = FALSE, success = FALSE;
76
77 /* Check if we have input parameters */
78 useropt = argc;
79 if(useropt < 2) {
80 printf("LUSOL v%d.%d.%d.%d - Linear equation solver - Development edition.\n",
81 LUSOL_VERMAJOR, LUSOL_VERMINOR, LUSOL_RELEASE, LUSOL_BUILD);
82 printf("Usage: lusol [-p <type>] [-t <tolerance>] [-m <memscalar>]\n");
83 printf(" [-s] <MatrixFile> [<RhsFile>]\n");
84 printf("Options: -p <type> <0..3> Selects pivot type.\n");
85 printf(" -t <tolerance> <1..100> Selects diagonal tolerance.\n");
86 printf(" -m <memscalar> <%d..100> Memory allocation scalar.\n", LUSOL_MULT_nz_a);
87 printf(" -b Solves using btran (rather than ftran).\n");
88 printf(" -r <times> Randomly replace a column and resolve.\n");
89 printf(" -s Show the computed solution.\n");
90 printf(" -blas <lib> Activates external optimized BLAS library.\n");
91 printf("Formats: Conventional RCV .TXT, MatrixMarket .MTX, or Harwell-Boeing .RUA\n");
92 printf("Author: Michael A. Saunders (original Fortran version)\n");
93 printf(" Kjell Eikland (modified C version)\n");
94 return;
95 }
96
97 /* Create the LUSOL object and set user options */
98 LUSOL = LUSOL_create(outunit, 0, LUSOL_PIVMOD_TPP, 0);
99 #if 1
100 LUSOL->luparm[LUSOL_IP_ACCELERATION] = LUSOL_OTHERORDER | LUSOL_ACCELERATE_L0;
101 #elif 0
102 LUSOL->luparm[LUSOL_IP_ACCELERATION] = LUSOL_AUTOORDER | LUSOL_ACCELERATE_L0;
103 #endif
104 LUSOL->luparm[LUSOL_IP_SCALAR_NZA] = 10;
105 i = 1;
106 n = 0;
107 filename[0] = '\0';
108 blasname[0] = '\0';
109 while((n == 0) && (i < argc)) {
110 if(strcmp("-p", argv[i]) == 0) {
111 i1 = i+1;
112 if((i1 < argc) && isNum(argv[i1][0])) {
113 i = i1;
114 m = atoi(argv[i]);
115 if(m < 0 || m > LUSOL_PIVMOD_MAX)
116 continue;
117 LUSOL->luparm[LUSOL_IP_PIVOTTYPE] = m;
118 }
119 }
120 else if(strcmp("-t", argv[i]) == 0) {
121 i1 = i+1;
122 if((i1 < argc) && isNum(argv[i1][0])) {
123 i = i1;
124 Amax = atof(argv[i]);
125 if(Amax < 1 || Amax > 100)
126 continue;
127 LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij] = Amax;
128 }
129 }
130 else if(strcmp("-m", argv[i]) == 0) {
131 i1 = i+1;
132 if((i1 < argc) && isNum(argv[i1][0])) {
133 i = i1;
134 m = atoi(argv[i]);
135 if(m < LUSOL_MULT_nz_a || m > 100)
136 continue;
137 LUSOL->luparm[LUSOL_IP_SCALAR_NZA] = atoi(argv[i]);
138 }
139 }
140 else if(strcmp("-s", argv[i]) == 0)
141 printsolution = TRUE;
142 else if(strcmp("-b", argv[i]) == 0)
143 ftran = FALSE;
144 else if(strcmp("-r", argv[i]) == 0) {
145 i1 = i+1;
146 if((i1 < argc) && isNum(argv[i1][0])) {
147 i = i1;
148 m = atoi(argv[i]);
149 if(m < 0 || m > 10)
150 continue;
151 }
152 else
153 m = 1;
154 srand((unsigned) time( NULL ));
155 replace = 2*m;
156 }
157 else if(strcmp("-blas", argv[i]) == 0) {
158 i1 = i+1;
159 if((i1 < argc) && getFileName(blasname, argv[i1])) {
160 if(!load_BLAS(blasname))
161 fprintf(outunit, "Could not load external BLAS library '%s'\n", blasname);
162 i = i1;
163 }
164 else
165 fprintf(outunit, "Ignoring incomplete parameter %d '%s'\n", i, argv[i]);
166 }
167 else {
168 if(getFileName(filename, argv[i])) {
169 useropt = i;
170 break;
171 }
172 else
173 fprintf(outunit, "Ignoring unknown parameter %d '%s'\n", i, argv[i]);
174 }
175 i++;
176 }
177
178 /* Obtain file extension and see if we must estimate matrix data size */
179 strcpy(fileext, strchr(argv[useropt], '.'));
180 /* Yin Zhang noted that _strupr() is not available on many Unix platforms. */
181 _strupr_(fileext);
182 /* _strupr(fileext);*/
183
184 /* Read conventional text file format */
185 if(strcmp(fileext, ".TXT") == 0) {
186 if(!ctf_size_A(filename, &maxm, &maxn, &maxnz))
187 goto x900;
188 }
189 /* Read MatrixMarket file format */
190 else if(strcmp(fileext, ".MTX") == 0) {
191 if(!mmf_size_A(filename, &maxm, &maxn, &maxnz))
192 goto x900;
193 }
194 /* Read Harwell-Boeing file format */
195 else if(strcmp(fileext, ".RUA") == 0) {
196 if(!hbf_size_A(filename, &maxm, &maxn, &maxnz))
197 goto x900;
198 }
199 else {
200 fprintf(outunit, "Unrecognized matrix file extension %s\n", fileext);
201 goto x900;
202 }
203
204 /* Create the arrays */
205
206 Aij = (REAL *) LUSOL_CALLOC(maxnz + BLAS_BASE, sizeof(REAL));
207 iA = (int *) LUSOL_CALLOC(maxnz + BLAS_BASE, sizeof(int));
208 jA = (int *) LUSOL_CALLOC(maxnz + BLAS_BASE, sizeof(int));
209 if(ftran)
210 lenb = maxm;
211 else
212 lenb = maxn;
213 b = (REAL *) LUSOL_CALLOC(lenb+BLAS_BASE, sizeof(REAL));
214 rhs = (REAL *) LUSOL_CALLOC(lenb+BLAS_BASE, sizeof(REAL));
215
216 if(ftran)
217 lenx = maxn;
218 else
219 lenx = maxm;
220 xexact = (REAL *) LUSOL_CALLOC(lenx+BLAS_BASE, sizeof(REAL));
221 r = (REAL *) LUSOL_CALLOC(lenx+BLAS_BASE, sizeof(REAL));
222 x = (REAL *) LUSOL_CALLOC(lenx+BLAS_BASE, sizeof(REAL));
223
224 /* -----------------------------------------------------------------
225 Read data files.
226 ----------------------------------------------------------------- */
227 fprintf(stdout, "\n========================================\n");
228 fprintf(stdout, "LUSOL v%d.%d.%d.%d - Linear equation solver",
229 LUSOL_VERMAJOR, LUSOL_VERMINOR, LUSOL_RELEASE, LUSOL_BUILD);
230 fprintf(stdout, "\n========================================\n");
231
232 /* -----------------------------------------------------------------
233 Read data for A
234 ----------------------------------------------------------------- */
235 /* Read conventional text file format */
236 if(strcmp(fileext, ".TXT") == 0) {
237 if(!ctf_read_A(argv[useropt], maxm, maxn, maxnz,
238 &m, &n, &nnzero, iA, jA, Aij))
239 goto x900;
240 }
241 /* Read MatrixMarket file format */
242 else if(strcmp(fileext, ".MTX") == 0) {
243 if(!mmf_read_A(argv[useropt], maxm, maxn, maxnz,
244 &m, &n, &nnzero, iA, jA, Aij))
245 goto x900;
246 }
247 /* Read Harwell-Boeing file format */
248 else if(strcmp(fileext, ".RUA") == 0) {
249 if(!hbf_read_A(argv[useropt], maxm, maxn, maxnz,
250 &m, &n, &nnzero, iA, jA, Aij))
251 goto x900;
252 }
253 else {
254 fprintf(outunit, "Error: Unrecognized matrix file extension %s\n", fileext);
255 goto x900;
256 }
257
258 /* -----------------------------------------------------------------
259 Read data for b
260 ----------------------------------------------------------------- */
261 /* Handle Harwell-Boeing case where the file contains a RHS */
262 i = strcmp(fileext, ".RUA");
263
264 if((useropt == argc-1) && (i != 0)) {
265 srand(timeNow());
266 i1 = m;
267 while(i1 > 0) {
268 test = rand();
269 i = RAND_MAX;
270 i = (int) ((test/i)*(m-1));
271 /* b[i+1] = 1.0; */
272 b[i+1] = i - 5;
273 i1--;
274 }
275 if(printsolution)
276 blockWriteREAL(outunit, "\nGenerated RHS vector", b, 1, lenb);
277 }
278 else {
279 if(i != 0)
280 useropt++;
281 strcpy(fileext, strchr(argv[useropt], '.'));
282 _strupr_(fileext);
283
284 /* Read conventional text file format */
285 if(strcmp(fileext, ".TXT") == 0) {
286 if(!ctf_read_b(argv[useropt], lenb, b))
287 goto x900;
288 }
289 /* Read MatrixMarket file format */
290 else if(strcmp(fileext, ".MTX") == 0) {
291 if(!mmf_read_b(argv[useropt], lenb, b))
292 goto x900;
293 }
294 /* Read Harwell-Boeing file format */
295 else if(strcmp(fileext, ".RUA") == 0) {
296 if(!hbf_read_b(argv[useropt], lenb, b))
297 goto x900;
298 }
299 else {
300 fprintf(outunit, "Error: Unrecognized vector file extension %s\n", fileext);
301 goto x900;
302 }
303 }
304 success = TRUE;
305
306 /* -----------------------------------------------------------------
307 Show data on input
308 ----------------------------------------------------------------- */
309 fprintf(outunit, "\nData read from:\n%s\n", filename);
310 test = (double) nnzero / ((double) m * (double) n);
311 test *= 100.0;
312 fprintf(outunit, "Rows = %d Columns = %d Non-zeros = %d Density =%8.4f%%\n",
313 m, n, nnzero, test);
314
315 /* -----------------------------------------------------------------
316 Load A into (a, indc, indr).
317 ----------------------------------------------------------------- */
318 #if 0 /* BUG !!! */
319 if(n != m)
320 LUSOL->luparm[LUSOL_IP_KEEPLU] = FALSE;
321 #endif
322 #ifdef LegacyTesting
323 LUSOL->luparm[LUSOL_IP_SCALAR_NZA] = LUSOL_MULT_nz_a;
324 LUSOL_sizeto(LUSOL, MAXROWS, MAXCOLS, MAXNZ*LUSOL_MULT_nz_a);
325 #endif
326
327 if(!LUSOL_assign(LUSOL, iA, jA, Aij, nnzero, TRUE)) {
328 fprintf(outunit, "Error: LUSOL failed due to insufficient memory.\n");
329 goto x900;
330 }
331
332 /* ------------------------------------------------------------------
333 Factor A = L U.
334 ------------------------------------------------------------------ */
335 nelem = nnzero;
336 inform = LUSOL_factorize( LUSOL);
337 if(inform > LUSOL_INFORM_SERIOUS) {
338 fprintf(outunit, "Error:\n%s\n", LUSOL_informstr(LUSOL, inform));
339 goto x900;
340 }
341 if(n != m)
342 goto x800;
343
344 /* Get the largest element in A; we use it below as an estimate
345 of ||A||_inf, even though it isn't a proper norm. */
346 Amax = LUSOL->parmlu[LUSOL_RP_MAXELEM_A];
347
348 /* ------------------------------------------------------------------
349 SOLVE A x = b or x'A = b'.
350 Save b first because lu6sol() overwrites the rhs.
351 ------------------------------------------------------------------ */
352 MEMCOPY(x, b, lenb+BLAS_BASE);
353
354 Resolve:
355 if(ftran)
356 inform = LUSOL_ftran(LUSOL, x, NULL, FALSE);
357 else
358 inform = LUSOL_btran(LUSOL, x, NULL);
359 if(inform > LUSOL_INFORM_SERIOUS) {
360 fprintf(outunit, "Error:\n%s\n", LUSOL_informstr(LUSOL, inform));
361 goto x900;
362 }
363 if(printsolution)
364 blockWriteREAL(outunit, "\nSolution vector", x, 1, lenb);
365
366 /* ------------------------------------------------------------------
367 Set r = b - Ax.
368 Find norm of r and x.
369 ------------------------------------------------------------------ */
370 MEMCOPY(r, b, lenb+BLAS_BASE);
371 for(k = 1; k <= nnzero; k++) {
372 i = iA[k]; /* Row number */
373 j = jA[k]; /* Column number */
374 if(ftran)
375 r[i] -= Aij[k]*x[j];
376 else
377 r[j] -= Aij[k]*x[i];
378 }
379 /* blockWriteREAL(outunit, "\nResidual vector", r, 1, lenb);*/
380 bnorm = dnormi( lenb, b );
381 rnorm = dnormi( lenb, r );
382 xnorm = dnormi( lenx, x );
383
384 /* ------------------------------------------------------------------
385 Report the findings.
386 ------------------------------------------------------------------ */
387 if(randcol > 0)
388 fprintf(outunit, "\n\nColumn %d was %s\n",
389 randcol,
390 (mod(replace,2) == 1 ? "replaced with random data" : "restored"));
391
392 x800:
393 fprintf(outunit, "\nLU size statistics (%d reallocations):\n",
394 LUSOL->expanded_a);
395 test = LUSOL->luparm[LUSOL_IP_NONZEROS_U0]+LUSOL->luparm[LUSOL_IP_NONZEROS_L0];
396 fprintf(outunit, "L0-size = %d U0-size = %d LU-nonzeros = %d Fill-in = %.1fx\n",
397 LUSOL->luparm[LUSOL_IP_NONZEROS_L0],
398 LUSOL->luparm[LUSOL_IP_NONZEROS_U0],
399 (int) test, test/nnzero);
400 if(n != m) {
401 fprintf(outunit, "%s with a factor tol. of %g identified %d singularities.\n",
402 LUSOL_pivotLabel(LUSOL), LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij],
403 LUSOL->luparm[LUSOL_IP_SINGULARITIES]);
404 goto x900;
405 }
406
407 test = rnorm / (Amax*xnorm);
408 fprintf(outunit, "\nAccuracy statistics:\n");
409 fprintf(outunit, "%s with a factor tol. of %g gave a rel.error of %g and %d singularities.\n",
410 LUSOL_pivotLabel(LUSOL), LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij], test,
411 LUSOL->luparm[LUSOL_IP_SINGULARITIES]);
412 fprintf(outunit, "Amax = %g bnorm = %g rnorm = %g xnorm = %g\n",
413 Amax, bnorm, rnorm, xnorm);
414
415 fprintf(outunit, "\n");
416
417 if (test <= 1.0e-8)
418 fprintf(outunit, "The equations were solved with very high accuracy.\n");
419 else if (test <= 1.0e-6)
420 fprintf(outunit, "The equations were solved with reasonably good accuracy.\n");
421 else {
422 if (test <= 1.0e-4)
423 fprintf(outunit, "Questionable accuracy; the LU factors may not be good enough.\n");
424 else
425 fprintf(outunit, "Poor accuracy; the LU factorization probably failed.\n");
426 if(LUSOL->luparm[LUSOL_IP_PIVOTTYPE] == LUSOL_PIVMOD_TRP)
427 fprintf(outunit, "Try a smaller factor tolerance (current is %g).\n",
428 LUSOL->parmlu[LUSOL_RP_FACTORMAX_Lij]);
429 else
430 fprintf(outunit, "Try a smaller factor tolerance and/or TRP pivoting.\n");
431 }
432
433 /* Check if we should replace a column and resolve */
434 if(replace > 0) {
435 replace--;
436 MEMCLEAR(x, lenb+BLAS_BASE);
437 if(mod(replace, 2) == 1) {
438 /* Randomly find a column and replace the column values with random data */
439 rnorm = rand();
440 randcol = (int) (n * rnorm / (RAND_MAX+1.0)) + 1;
441 for(i = 1; i < m; i++)
442 x[i] = Amax * rand() / RAND_MAX;
443 }
444 else {
445 /* Put the previously replaced column back and resolve */
446 for (k = 1; k <= nnzero; k++) {
447 i = iA[k];
448 j = jA[k];
449 if(j == randcol)
450 x[i] = Aij[k];
451 }
452 }
453 inform = LUSOL_replaceColumn(LUSOL, randcol, x);
454 MEMCOPY(b, x, lenb+BLAS_BASE);
455 if(inform != LUSOL_INFORM_LUSUCCESS)
456 fprintf(outunit, "Error:\n%s\n", LUSOL_informstr(LUSOL, inform));
457 else
458 goto Resolve;
459 }
460
461
462 /* Free memory */
463 x900:
464 if(!success)
465 fprintf(outunit, "Insufficient memory or data file not found.\n");
466
467 LUSOL_FREE(Aij);
468 LUSOL_FREE(b);
469 LUSOL_FREE(xexact);
470 LUSOL_FREE(iA);
471 LUSOL_FREE(jA);
472
473 LUSOL_FREE(rhs);
474 LUSOL_FREE(r);
475 LUSOL_FREE(x);
476
477 #if 0
478 LUSOL_dump(NULL, LUSOL);
479 -blas "atlas_AXP_512_360.dll" -b -p 1 "STP3D_A.MTX"
480 "C:\Shared Files\Visual Studio Projects\LU\MatrixMarket\sherman5.mtx" "C:\Shared Files\Visual Studio Projects\LU\MatrixMarket\sherman5_rhs1.mtx"
481 A6805.txt b6805.txt
482 A10009.txt b10009.txt
483 fidap005.mtx fidap005_rhs1.mtx
484 fidapm05.mtx fidapm05_rhs1.mtx
485 -b -p 1 "basis.mtx"
486 -b -p 1 "LU-test3.mtx"
487 #endif
488
489 LUSOL_free(LUSOL);
490
491 /* End of main program for Test of LUSOL */
492 }
493
494