1 /*! \file
2 Copyright (c) 2003, The Regents of the University of California, through
3 Lawrence Berkeley National Laboratory (subject to receipt of any required
4 approvals from U.S. Dept. of Energy)
5
6 All rights reserved.
7
8 The source code is distributed under BSD license, see the file License.txt
9 at the top-level directory.
10 */
11
12 /*! @file dutil.c
13 * \brief Matrix utility functions
14 *
15 * <pre>
16 * -- SuperLU routine (version 3.1) --
17 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
18 * and Lawrence Berkeley National Lab.
19 * August 1, 2008
20 *
21 * Copyright (c) 1994 by Xerox Corporation. All rights reserved.
22 *
23 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
24 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
25 *
26 * Permission is hereby granted to use or copy this program for any
27 * purpose, provided the above notices are retained on all copies.
28 * Permission to modify the code and to distribute modified code is
29 * granted, provided the above notices are retained, and a notice that
30 * the code was modified is included with the above copyright notice.
31 * </pre>
32 */
33
34
35 #include <math.h>
36 #include "slu_ddefs.h"
37
38 void
dCreate_CompCol_Matrix(SuperMatrix * A,int m,int n,int nnz,double * nzval,int * rowind,int * colptr,Stype_t stype,Dtype_t dtype,Mtype_t mtype)39 dCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int nnz,
40 double *nzval, int *rowind, int *colptr,
41 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
42 {
43 NCformat *Astore;
44
45 A->Stype = stype;
46 A->Dtype = dtype;
47 A->Mtype = mtype;
48 A->nrow = m;
49 A->ncol = n;
50 A->Store = (void *) SUPERLU_MALLOC( sizeof(NCformat) );
51 if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
52 Astore = A->Store;
53 Astore->nnz = nnz;
54 Astore->nzval = nzval;
55 Astore->rowind = rowind;
56 Astore->colptr = colptr;
57 }
58
59 void
dCreate_CompRow_Matrix(SuperMatrix * A,int m,int n,int nnz,double * nzval,int * colind,int * rowptr,Stype_t stype,Dtype_t dtype,Mtype_t mtype)60 dCreate_CompRow_Matrix(SuperMatrix *A, int m, int n, int nnz,
61 double *nzval, int *colind, int *rowptr,
62 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
63 {
64 NRformat *Astore;
65
66 A->Stype = stype;
67 A->Dtype = dtype;
68 A->Mtype = mtype;
69 A->nrow = m;
70 A->ncol = n;
71 A->Store = (void *) SUPERLU_MALLOC( sizeof(NRformat) );
72 if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
73 Astore = A->Store;
74 Astore->nnz = nnz;
75 Astore->nzval = nzval;
76 Astore->colind = colind;
77 Astore->rowptr = rowptr;
78 }
79
80 /*! \brief Copy matrix A into matrix B. */
81 void
dCopy_CompCol_Matrix(SuperMatrix * A,SuperMatrix * B)82 dCopy_CompCol_Matrix(SuperMatrix *A, SuperMatrix *B)
83 {
84 NCformat *Astore, *Bstore;
85 int ncol, nnz, i;
86
87 B->Stype = A->Stype;
88 B->Dtype = A->Dtype;
89 B->Mtype = A->Mtype;
90 B->nrow = A->nrow;;
91 B->ncol = ncol = A->ncol;
92 Astore = (NCformat *) A->Store;
93 Bstore = (NCformat *) B->Store;
94 Bstore->nnz = nnz = Astore->nnz;
95 for (i = 0; i < nnz; ++i)
96 ((double *)Bstore->nzval)[i] = ((double *)Astore->nzval)[i];
97 for (i = 0; i < nnz; ++i) Bstore->rowind[i] = Astore->rowind[i];
98 for (i = 0; i <= ncol; ++i) Bstore->colptr[i] = Astore->colptr[i];
99 }
100
101
102 void
dCreate_Dense_Matrix(SuperMatrix * X,int m,int n,double * x,int ldx,Stype_t stype,Dtype_t dtype,Mtype_t mtype)103 dCreate_Dense_Matrix(SuperMatrix *X, int m, int n, double *x, int ldx,
104 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
105 {
106 DNformat *Xstore;
107
108 X->Stype = stype;
109 X->Dtype = dtype;
110 X->Mtype = mtype;
111 X->nrow = m;
112 X->ncol = n;
113 X->Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
114 if ( !(X->Store) ) ABORT("SUPERLU_MALLOC fails for X->Store");
115 Xstore = (DNformat *) X->Store;
116 Xstore->lda = ldx;
117 Xstore->nzval = (double *) x;
118 }
119
120 void
dCopy_Dense_Matrix(int M,int N,double * X,int ldx,double * Y,int ldy)121 dCopy_Dense_Matrix(int M, int N, double *X, int ldx,
122 double *Y, int ldy)
123 {
124 /*! \brief Copies a two-dimensional matrix X to another matrix Y.
125 */
126 int i, j;
127
128 for (j = 0; j < N; ++j)
129 for (i = 0; i < M; ++i)
130 Y[i + j*ldy] = X[i + j*ldx];
131 }
132
133 void
dCreate_SuperNode_Matrix(SuperMatrix * L,int m,int n,int nnz,double * nzval,int * nzval_colptr,int * rowind,int * rowind_colptr,int * col_to_sup,int * sup_to_col,Stype_t stype,Dtype_t dtype,Mtype_t mtype)134 dCreate_SuperNode_Matrix(SuperMatrix *L, int m, int n, int nnz,
135 double *nzval, int *nzval_colptr, int *rowind,
136 int *rowind_colptr, int *col_to_sup, int *sup_to_col,
137 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
138 {
139 SCformat *Lstore;
140
141 L->Stype = stype;
142 L->Dtype = dtype;
143 L->Mtype = mtype;
144 L->nrow = m;
145 L->ncol = n;
146 L->Store = (void *) SUPERLU_MALLOC( sizeof(SCformat) );
147 if ( !(L->Store) ) ABORT("SUPERLU_MALLOC fails for L->Store");
148 Lstore = L->Store;
149 Lstore->nnz = nnz;
150 Lstore->nsuper = col_to_sup[n];
151 Lstore->nzval = nzval;
152 Lstore->nzval_colptr = nzval_colptr;
153 Lstore->rowind = rowind;
154 Lstore->rowind_colptr = rowind_colptr;
155 Lstore->col_to_sup = col_to_sup;
156 Lstore->sup_to_col = sup_to_col;
157
158 }
159
160
161 /*! \brief Convert a row compressed storage into a column compressed storage.
162 */
163 void
dCompRow_to_CompCol(int m,int n,int nnz,double * a,int * colind,int * rowptr,double ** at,int ** rowind,int ** colptr)164 dCompRow_to_CompCol(int m, int n, int nnz,
165 double *a, int *colind, int *rowptr,
166 double **at, int **rowind, int **colptr)
167 {
168 register int i, j, col, relpos;
169 int *marker;
170
171 /* Allocate storage for another copy of the matrix. */
172 *at = (double *) doubleMalloc(nnz);
173 *rowind = (int *) intMalloc(nnz);
174 *colptr = (int *) intMalloc(n+1);
175 marker = (int *) intCalloc(n);
176
177 /* Get counts of each column of A, and set up column pointers */
178 for (i = 0; i < m; ++i)
179 for (j = rowptr[i]; j < rowptr[i+1]; ++j) ++marker[colind[j]];
180 (*colptr)[0] = 0;
181 for (j = 0; j < n; ++j) {
182 (*colptr)[j+1] = (*colptr)[j] + marker[j];
183 marker[j] = (*colptr)[j];
184 }
185
186 /* Transfer the matrix into the compressed column storage. */
187 for (i = 0; i < m; ++i) {
188 for (j = rowptr[i]; j < rowptr[i+1]; ++j) {
189 col = colind[j];
190 relpos = marker[col];
191 (*rowind)[relpos] = i;
192 (*at)[relpos] = a[j];
193 ++marker[col];
194 }
195 }
196
197 SUPERLU_FREE(marker);
198 }
199
200
201 void
dPrint_CompCol_Matrix(char * what,SuperMatrix * A)202 dPrint_CompCol_Matrix(char *what, SuperMatrix *A)
203 {
204 NCformat *Astore;
205 register int i,n;
206 double *dp;
207
208 printf("\nCompCol matrix %s:\n", what);
209 printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
210 n = A->ncol;
211 Astore = (NCformat *) A->Store;
212 dp = (double *) Astore->nzval;
213 printf("nrow %d, ncol %d, nnz %d\n", A->nrow,A->ncol,Astore->nnz);
214 printf("nzval: ");
215 for (i = 0; i < Astore->colptr[n]; ++i) printf("%f ", dp[i]);
216 printf("\nrowind: ");
217 for (i = 0; i < Astore->colptr[n]; ++i) printf("%d ", Astore->rowind[i]);
218 printf("\ncolptr: ");
219 for (i = 0; i <= n; ++i) printf("%d ", Astore->colptr[i]);
220 printf("\n");
221 fflush(stdout);
222 }
223
224 void
dPrint_SuperNode_Matrix(char * what,SuperMatrix * A)225 dPrint_SuperNode_Matrix(char *what, SuperMatrix *A)
226 {
227 SCformat *Astore;
228 register int i, j, k, c, d, n, nsup;
229 double *dp;
230 int *col_to_sup, *sup_to_col, *rowind, *rowind_colptr;
231
232 printf("\nSuperNode matrix %s:\n", what);
233 printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
234 n = A->ncol;
235 Astore = (SCformat *) A->Store;
236 dp = (double *) Astore->nzval;
237 col_to_sup = Astore->col_to_sup;
238 sup_to_col = Astore->sup_to_col;
239 rowind_colptr = Astore->rowind_colptr;
240 rowind = Astore->rowind;
241 printf("nrow %d, ncol %d, nnz %d, nsuper %d\n",
242 A->nrow,A->ncol,Astore->nnz,Astore->nsuper);
243 printf("nzval:\n");
244 for (k = 0; k <= Astore->nsuper; ++k) {
245 c = sup_to_col[k];
246 nsup = sup_to_col[k+1] - c;
247 for (j = c; j < c + nsup; ++j) {
248 d = Astore->nzval_colptr[j];
249 for (i = rowind_colptr[c]; i < rowind_colptr[c+1]; ++i) {
250 printf("%d\t%d\t%e\n", rowind[i], j, dp[d++]);
251 }
252 }
253 }
254 #if 0
255 for (i = 0; i < Astore->nzval_colptr[n]; ++i) printf("%f ", dp[i]);
256 #endif
257 printf("\nnzval_colptr: ");
258 for (i = 0; i <= n; ++i) printf("%d ", Astore->nzval_colptr[i]);
259 printf("\nrowind: ");
260 for (i = 0; i < Astore->rowind_colptr[n]; ++i)
261 printf("%d ", Astore->rowind[i]);
262 printf("\nrowind_colptr: ");
263 for (i = 0; i <= n; ++i) printf("%d ", Astore->rowind_colptr[i]);
264 printf("\ncol_to_sup: ");
265 for (i = 0; i < n; ++i) printf("%d ", col_to_sup[i]);
266 printf("\nsup_to_col: ");
267 for (i = 0; i <= Astore->nsuper+1; ++i)
268 printf("%d ", sup_to_col[i]);
269 printf("\n");
270 fflush(stdout);
271 }
272
273 void
dPrint_Dense_Matrix(char * what,SuperMatrix * A)274 dPrint_Dense_Matrix(char *what, SuperMatrix *A)
275 {
276 DNformat *Astore = (DNformat *) A->Store;
277 register int i, j, lda = Astore->lda;
278 double *dp;
279
280 printf("\nDense matrix %s:\n", what);
281 printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
282 dp = (double *) Astore->nzval;
283 printf("nrow %d, ncol %d, lda %d\n", A->nrow,A->ncol,lda);
284 printf("\nnzval: ");
285 for (j = 0; j < A->ncol; ++j) {
286 for (i = 0; i < A->nrow; ++i) printf("%f ", dp[i + j*lda]);
287 printf("\n");
288 }
289 printf("\n");
290 fflush(stdout);
291 }
292
293 /*! \brief Diagnostic print of column "jcol" in the U/L factor.
294 */
295 void
dprint_lu_col(char * msg,int jcol,int pivrow,int * xprune,GlobalLU_t * Glu)296 dprint_lu_col(char *msg, int jcol, int pivrow, int *xprune, GlobalLU_t *Glu)
297 {
298 int i, k, fsupc;
299 int *xsup, *supno;
300 int *xlsub, *lsub;
301 double *lusup;
302 int *xlusup;
303 double *ucol;
304 int *usub, *xusub;
305
306 xsup = Glu->xsup;
307 supno = Glu->supno;
308 lsub = Glu->lsub;
309 xlsub = Glu->xlsub;
310 lusup = (double *) Glu->lusup;
311 xlusup = Glu->xlusup;
312 ucol = (double *) Glu->ucol;
313 usub = Glu->usub;
314 xusub = Glu->xusub;
315
316 printf("%s", msg);
317 printf("col %d: pivrow %d, supno %d, xprune %d\n",
318 jcol, pivrow, supno[jcol], xprune[jcol]);
319
320 printf("\tU-col:\n");
321 for (i = xusub[jcol]; i < xusub[jcol+1]; i++)
322 printf("\t%d%10.4f\n", usub[i], ucol[i]);
323 printf("\tL-col in rectangular snode:\n");
324 fsupc = xsup[supno[jcol]]; /* first col of the snode */
325 i = xlsub[fsupc];
326 k = xlusup[jcol];
327 while ( i < xlsub[fsupc+1] && k < xlusup[jcol+1] ) {
328 printf("\t%d\t%10.4f\n", lsub[i], lusup[k]);
329 i++; k++;
330 }
331 fflush(stdout);
332 }
333
334
335 /*! \brief Check whether tempv[] == 0. This should be true before and after calling any numeric routines, i.e., "panel_bmod" and "column_bmod".
336 */
dcheck_tempv(int n,double * tempv)337 void dcheck_tempv(int n, double *tempv)
338 {
339 int i;
340
341 for (i = 0; i < n; i++) {
342 if (tempv[i] != 0.0)
343 {
344 fprintf(stderr,"tempv[%d] = %f\n", i,tempv[i]);
345 ABORT("dcheck_tempv");
346 }
347 }
348 }
349
350
351 void
dGenXtrue(int n,int nrhs,double * x,int ldx)352 dGenXtrue(int n, int nrhs, double *x, int ldx)
353 {
354 int i, j;
355 for (j = 0; j < nrhs; ++j)
356 for (i = 0; i < n; ++i) {
357 x[i + j*ldx] = 1.0;/* + (double)(i+1.)/n;*/
358 }
359 }
360
361 /*! \brief Let rhs[i] = sum of i-th row of A, so the solution vector is all 1's
362 */
363 void
dFillRHS(trans_t trans,int nrhs,double * x,int ldx,SuperMatrix * A,SuperMatrix * B)364 dFillRHS(trans_t trans, int nrhs, double *x, int ldx,
365 SuperMatrix *A, SuperMatrix *B)
366 {
367 NCformat *Astore;
368 double *Aval;
369 DNformat *Bstore;
370 double *rhs;
371 double one = 1.0;
372 double zero = 0.0;
373 int ldc;
374 char transc[1];
375
376 Astore = A->Store;
377 Aval = (double *) Astore->nzval;
378 Bstore = B->Store;
379 rhs = Bstore->nzval;
380 ldc = Bstore->lda;
381
382 if ( trans == NOTRANS ) *(unsigned char *)transc = 'N';
383 else *(unsigned char *)transc = 'T';
384
385 sp_dgemm(transc, "N", A->nrow, nrhs, A->ncol, one, A,
386 x, ldx, zero, rhs, ldc);
387
388 }
389
390 /*! \brief Fills a double precision array with a given value.
391 */
392 void
dfill(double * a,int alen,double dval)393 dfill(double *a, int alen, double dval)
394 {
395 register int i;
396 for (i = 0; i < alen; i++) a[i] = dval;
397 }
398
399
400
401 /*! \brief Check the inf-norm of the error vector
402 */
dinf_norm_error(int nrhs,SuperMatrix * X,double * xtrue)403 void dinf_norm_error(int nrhs, SuperMatrix *X, double *xtrue)
404 {
405 DNformat *Xstore;
406 double err, xnorm;
407 double *Xmat, *soln_work;
408 int i, j;
409
410 Xstore = X->Store;
411 Xmat = Xstore->nzval;
412
413 for (j = 0; j < nrhs; j++) {
414 soln_work = &Xmat[j*Xstore->lda];
415 err = xnorm = 0.0;
416 for (i = 0; i < X->nrow; i++) {
417 err = SUPERLU_MAX(err, fabs(soln_work[i] - xtrue[i]));
418 xnorm = SUPERLU_MAX(xnorm, fabs(soln_work[i]));
419 }
420 err = err / xnorm;
421 printf("||X - Xtrue||/||X|| = %e\n", err);
422 }
423 }
424
425
426
427 /*! \brief Print performance of the code. */
428 void
dPrintPerf(SuperMatrix * L,SuperMatrix * U,mem_usage_t * mem_usage,double rpg,double rcond,double * ferr,double * berr,char * equed,SuperLUStat_t * stat)429 dPrintPerf(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage,
430 double rpg, double rcond, double *ferr,
431 double *berr, char *equed, SuperLUStat_t *stat)
432 {
433 SCformat *Lstore;
434 NCformat *Ustore;
435 double *utime;
436 flops_t *ops;
437
438 utime = stat->utime;
439 ops = stat->ops;
440
441 if ( utime[FACT] != 0. )
442 printf("Factor flops = %e\tMflops = %8.2f\n", ops[FACT],
443 ops[FACT]*1e-6/utime[FACT]);
444 printf("Identify relaxed snodes = %8.2f\n", utime[RELAX]);
445 if ( utime[SOLVE] != 0. )
446 printf("Solve flops = %.0f, Mflops = %8.2f\n", ops[SOLVE],
447 ops[SOLVE]*1e-6/utime[SOLVE]);
448
449 Lstore = (SCformat *) L->Store;
450 Ustore = (NCformat *) U->Store;
451 printf("\tNo of nonzeros in factor L = %d\n", Lstore->nnz);
452 printf("\tNo of nonzeros in factor U = %d\n", Ustore->nnz);
453 printf("\tNo of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz);
454
455 printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
456 mem_usage->for_lu/1e6, mem_usage->total_needed/1e6);
457 printf("Number of memory expansions: %d\n", stat->expansions);
458
459 printf("\tFactor\tMflops\tSolve\tMflops\tEtree\tEquil\tRcond\tRefine\n");
460 printf("PERF:%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f\n",
461 utime[FACT], ops[FACT]*1e-6/utime[FACT],
462 utime[SOLVE], ops[SOLVE]*1e-6/utime[SOLVE],
463 utime[ETREE], utime[EQUIL], utime[RCOND], utime[REFINE]);
464
465 printf("\tRpg\t\tRcond\t\tFerr\t\tBerr\t\tEquil?\n");
466 printf("NUM:\t%e\t%e\t%e\t%e\t%s\n",
467 rpg, rcond, ferr[0], berr[0], equed);
468
469 }
470
471
472
473
print_double_vec(char * what,int n,double * vec)474 print_double_vec(char *what, int n, double *vec)
475 {
476 int i;
477 printf("%s: n %d\n", what, n);
478 for (i = 0; i < n; ++i) printf("%d\t%f\n", i, vec[i]);
479 return 0;
480 }
481
482