Line data Source code
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
39 0 : dCreate_CompCol_Matrix(SuperMatrix *A, int m, int n, int_t nnz,
40 : double *nzval, int_t *rowind, int_t *colptr,
41 : Stype_t stype, Dtype_t dtype, Mtype_t mtype)
42 : {
43 : NCformat *Astore;
44 :
45 0 : A->Stype = stype;
46 0 : A->Dtype = dtype;
47 0 : A->Mtype = mtype;
48 0 : A->nrow = m;
49 0 : A->ncol = n;
50 0 : A->Store = (void *) SUPERLU_MALLOC( sizeof(NCformat) );
51 0 : if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
52 0 : Astore = A->Store;
53 0 : Astore->nnz = nnz;
54 0 : Astore->nzval = nzval;
55 0 : Astore->rowind = rowind;
56 0 : Astore->colptr = colptr;
57 0 : }
58 :
59 : void
60 0 : dCreate_CompRow_Matrix(SuperMatrix *A, int m, int n, int_t nnz,
61 : double *nzval, int_t *colind, int_t *rowptr,
62 : Stype_t stype, Dtype_t dtype, Mtype_t mtype)
63 : {
64 : NRformat *Astore;
65 :
66 0 : A->Stype = stype;
67 0 : A->Dtype = dtype;
68 0 : A->Mtype = mtype;
69 0 : A->nrow = m;
70 0 : A->ncol = n;
71 0 : A->Store = (void *) SUPERLU_MALLOC( sizeof(NRformat) );
72 0 : if ( !(A->Store) ) ABORT("SUPERLU_MALLOC fails for A->Store");
73 0 : Astore = A->Store;
74 0 : Astore->nnz = nnz;
75 0 : Astore->nzval = nzval;
76 0 : Astore->colind = colind;
77 0 : Astore->rowptr = rowptr;
78 0 : }
79 :
80 : /*! \brief Copy matrix A into matrix B. */
81 : void
82 0 : dCopy_CompCol_Matrix(SuperMatrix *A, SuperMatrix *B)
83 : {
84 : NCformat *Astore, *Bstore;
85 : int ncol, nnz, i;
86 :
87 0 : B->Stype = A->Stype;
88 0 : B->Dtype = A->Dtype;
89 0 : B->Mtype = A->Mtype;
90 0 : B->nrow = A->nrow;;
91 0 : B->ncol = ncol = A->ncol;
92 0 : Astore = (NCformat *) A->Store;
93 0 : Bstore = (NCformat *) B->Store;
94 0 : Bstore->nnz = nnz = Astore->nnz;
95 0 : for (i = 0; i < nnz; ++i)
96 0 : ((double *)Bstore->nzval)[i] = ((double *)Astore->nzval)[i];
97 0 : for (i = 0; i < nnz; ++i) Bstore->rowind[i] = Astore->rowind[i];
98 0 : for (i = 0; i <= ncol; ++i) Bstore->colptr[i] = Astore->colptr[i];
99 0 : }
100 :
101 :
102 : void
103 0 : 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 0 : X->Stype = stype;
109 0 : X->Dtype = dtype;
110 0 : X->Mtype = mtype;
111 0 : X->nrow = m;
112 0 : X->ncol = n;
113 0 : X->Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
114 0 : if ( !(X->Store) ) ABORT("SUPERLU_MALLOC fails for X->Store");
115 0 : Xstore = (DNformat *) X->Store;
116 0 : Xstore->lda = ldx;
117 0 : Xstore->nzval = (double *) x;
118 0 : }
119 :
120 : void
121 0 : 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 0 : for (j = 0; j < N; ++j)
129 0 : for (i = 0; i < M; ++i)
130 0 : Y[i + j*ldy] = X[i + j*ldx];
131 0 : }
132 :
133 : void
134 0 : dCreate_SuperNode_Matrix(SuperMatrix *L, int m, int n, int_t nnz,
135 : double *nzval, int_t *nzval_colptr, int_t *rowind,
136 : int_t *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 0 : L->Stype = stype;
142 0 : L->Dtype = dtype;
143 0 : L->Mtype = mtype;
144 0 : L->nrow = m;
145 0 : L->ncol = n;
146 0 : L->Store = (void *) SUPERLU_MALLOC( sizeof(SCformat) );
147 0 : if ( !(L->Store) ) ABORT("SUPERLU_MALLOC fails for L->Store");
148 0 : Lstore = L->Store;
149 0 : Lstore->nnz = nnz;
150 0 : Lstore->nsuper = col_to_sup[n];
151 0 : Lstore->nzval = nzval;
152 0 : Lstore->nzval_colptr = nzval_colptr;
153 0 : Lstore->rowind = rowind;
154 0 : Lstore->rowind_colptr = rowind_colptr;
155 0 : Lstore->col_to_sup = col_to_sup;
156 0 : Lstore->sup_to_col = sup_to_col;
157 :
158 0 : }
159 :
160 :
161 : /*! \brief Convert a row compressed storage into a column compressed storage.
162 : */
163 : void
164 0 : dCompRow_to_CompCol(int m, int n, int_t nnz,
165 : double *a, int_t *colind, int_t *rowptr,
166 : double **at, int_t **rowind, int_t **colptr)
167 : {
168 : register int i, j, col, relpos;
169 : int_t *marker;
170 :
171 : /* Allocate storage for another copy of the matrix. */
172 0 : *at = (double *) doubleMalloc(nnz);
173 0 : *rowind = (int_t *) intMalloc(nnz);
174 0 : *colptr = (int_t *) intMalloc(n+1);
175 0 : marker = (int_t *) intCalloc(n);
176 :
177 : /* Get counts of each column of A, and set up column pointers */
178 0 : for (i = 0; i < m; ++i)
179 0 : for (j = rowptr[i]; j < rowptr[i+1]; ++j) ++marker[colind[j]];
180 0 : (*colptr)[0] = 0;
181 0 : for (j = 0; j < n; ++j) {
182 0 : (*colptr)[j+1] = (*colptr)[j] + marker[j];
183 0 : marker[j] = (*colptr)[j];
184 : }
185 :
186 : /* Transfer the matrix into the compressed column storage. */
187 0 : for (i = 0; i < m; ++i) {
188 0 : for (j = rowptr[i]; j < rowptr[i+1]; ++j) {
189 0 : col = colind[j];
190 0 : relpos = marker[col];
191 0 : (*rowind)[relpos] = i;
192 0 : (*at)[relpos] = a[j];
193 0 : ++marker[col];
194 : }
195 : }
196 :
197 0 : SUPERLU_FREE(marker);
198 0 : }
199 :
200 :
201 : void
202 0 : dPrint_CompCol_Matrix(char *what, SuperMatrix *A)
203 : {
204 : NCformat *Astore;
205 : register int_t i;
206 : register int n;
207 : double *dp;
208 :
209 : printf("\nCompCol matrix %s:\n", what);
210 0 : printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
211 0 : n = A->ncol;
212 0 : Astore = (NCformat *) A->Store;
213 0 : dp = (double *) Astore->nzval;
214 0 : printf("nrow %d, ncol %d, nnz %ld\n", (int)A->nrow, (int)A->ncol, (long)Astore->nnz);
215 : printf("nzval: ");
216 0 : for (i = 0; i < Astore->colptr[n]; ++i) printf("%f ", dp[i]);
217 : printf("\nrowind: ");
218 0 : for (i = 0; i < Astore->colptr[n]; ++i) printf("%ld ", (long)Astore->rowind[i]);
219 : printf("\ncolptr: ");
220 0 : for (i = 0; i <= n; ++i) printf("%ld ", (long)Astore->colptr[i]);
221 : printf("\n");
222 0 : fflush(stdout);
223 0 : }
224 :
225 : void
226 0 : dPrint_SuperNode_Matrix(char *what, SuperMatrix *A)
227 : {
228 : SCformat *Astore;
229 : register int_t i, j, k, c, d, n, nsup;
230 : double *dp;
231 : int *col_to_sup, *sup_to_col;
232 : int_t *rowind, *rowind_colptr;
233 :
234 : printf("\nSuperNode matrix %s:\n", what);
235 0 : printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
236 0 : n = A->ncol;
237 0 : Astore = (SCformat *) A->Store;
238 0 : dp = (double *) Astore->nzval;
239 0 : col_to_sup = Astore->col_to_sup;
240 0 : sup_to_col = Astore->sup_to_col;
241 0 : rowind_colptr = Astore->rowind_colptr;
242 0 : rowind = Astore->rowind;
243 0 : printf("nrow %d, ncol %d, nnz %lld, nsuper %d\n",
244 0 : (int)A->nrow, (int)A->ncol, (long long) Astore->nnz, (int)Astore->nsuper);
245 : printf("nzval:\n");
246 0 : for (k = 0; k <= Astore->nsuper; ++k) {
247 0 : c = sup_to_col[k];
248 0 : nsup = sup_to_col[k+1] - c;
249 0 : for (j = c; j < c + nsup; ++j) {
250 0 : d = Astore->nzval_colptr[j];
251 0 : for (i = rowind_colptr[c]; i < rowind_colptr[c+1]; ++i) {
252 0 : printf("%d\t%d\t%e\n", (int)rowind[i], (int)j, dp[d++]);
253 : }
254 : }
255 : }
256 : #if 0
257 : for (i = 0; i < Astore->nzval_colptr[n]; ++i) printf("%f ", dp[i]);
258 : #endif
259 : printf("\nnzval_colptr: ");
260 0 : for (i = 0; i <= n; ++i) printf("%lld ", (long long)Astore->nzval_colptr[i]);
261 : printf("\nrowind: ");
262 0 : for (i = 0; i < Astore->rowind_colptr[n]; ++i)
263 0 : printf("%lld ", (long long)Astore->rowind[i]);
264 : printf("\nrowind_colptr: ");
265 0 : for (i = 0; i <= n; ++i) printf("%lld ", (long long)Astore->rowind_colptr[i]);
266 : printf("\ncol_to_sup: ");
267 0 : for (i = 0; i < n; ++i) printf("%d ", col_to_sup[i]);
268 : printf("\nsup_to_col: ");
269 0 : for (i = 0; i <= Astore->nsuper+1; ++i)
270 0 : printf("%d ", sup_to_col[i]);
271 : printf("\n");
272 0 : fflush(stdout);
273 0 : }
274 :
275 : void
276 0 : dPrint_Dense_Matrix(char *what, SuperMatrix *A)
277 : {
278 0 : DNformat *Astore = (DNformat *) A->Store;
279 0 : register int i, j, lda = Astore->lda;
280 : double *dp;
281 :
282 : printf("\nDense matrix %s:\n", what);
283 0 : printf("Stype %d, Dtype %d, Mtype %d\n", A->Stype,A->Dtype,A->Mtype);
284 0 : dp = (double *) Astore->nzval;
285 0 : printf("nrow %d, ncol %d, lda %d\n", (int)A->nrow, (int)A->ncol, lda);
286 : printf("\nnzval: ");
287 0 : for (j = 0; j < A->ncol; ++j) {
288 0 : for (i = 0; i < A->nrow; ++i) printf("%f ", dp[i + j*lda]);
289 : printf("\n");
290 : }
291 : printf("\n");
292 0 : fflush(stdout);
293 0 : }
294 :
295 : /*! \brief Diagnostic print of column "jcol" in the U/L factor.
296 : */
297 : void
298 0 : dprint_lu_col(char *msg, int jcol, int pivrow, int_t *xprune, GlobalLU_t *Glu)
299 : {
300 : int_t i, k;
301 : int *xsup, *supno, fsupc;
302 : int_t *xlsub, *lsub;
303 : double *lusup;
304 : int_t *xlusup;
305 : double *ucol;
306 : int_t *usub, *xusub;
307 :
308 0 : xsup = Glu->xsup;
309 0 : supno = Glu->supno;
310 0 : lsub = Glu->lsub;
311 0 : xlsub = Glu->xlsub;
312 0 : lusup = (double *) Glu->lusup;
313 0 : xlusup = Glu->xlusup;
314 0 : ucol = (double *) Glu->ucol;
315 0 : usub = Glu->usub;
316 0 : xusub = Glu->xusub;
317 :
318 : printf("%s", msg);
319 0 : printf("col %d: pivrow %d, supno %d, xprune %lld\n",
320 0 : jcol, pivrow, supno[jcol], (long long) xprune[jcol]);
321 :
322 : printf("\tU-col:\n");
323 0 : for (i = xusub[jcol]; i < xusub[jcol+1]; i++)
324 0 : printf("\t%d%10.4f\n", (int)usub[i], ucol[i]);
325 : printf("\tL-col in rectangular snode:\n");
326 0 : fsupc = xsup[supno[jcol]]; /* first col of the snode */
327 0 : i = xlsub[fsupc];
328 0 : k = xlusup[jcol];
329 0 : while ( i < xlsub[fsupc+1] && k < xlusup[jcol+1] ) {
330 0 : printf("\t%d\t%10.4f\n", (int)lsub[i], lusup[k]);
331 0 : i++; k++;
332 : }
333 0 : fflush(stdout);
334 0 : }
335 :
336 :
337 : /*! \brief Check whether tempv[] == 0. This should be true before and after calling any numeric routines, i.e., "panel_bmod" and "column_bmod".
338 : */
339 0 : void dcheck_tempv(int n, double *tempv)
340 : {
341 : int i;
342 :
343 0 : for (i = 0; i < n; i++) {
344 0 : if (tempv[i] != 0.0)
345 : {
346 0 : fprintf(stderr,"tempv[%d] = %f\n", i,tempv[i]);
347 0 : ABORT("dcheck_tempv");
348 : }
349 : }
350 0 : }
351 :
352 :
353 : void
354 0 : dGenXtrue(int n, int nrhs, double *x, int ldx)
355 : {
356 : int i, j;
357 0 : for (j = 0; j < nrhs; ++j)
358 0 : for (i = 0; i < n; ++i) {
359 0 : x[i + j*ldx] = 1.0;/* + (double)(i+1.)/n;*/
360 : }
361 0 : }
362 :
363 : /*! \brief Let rhs[i] = sum of i-th row of A, so the solution vector is all 1's
364 : */
365 : void
366 0 : dFillRHS(trans_t trans, int nrhs, double *x, int ldx,
367 : SuperMatrix *A, SuperMatrix *B)
368 : {
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 0 : Bstore = B->Store;
379 0 : rhs = Bstore->nzval;
380 0 : ldc = Bstore->lda;
381 :
382 0 : if ( trans == NOTRANS ) *(unsigned char *)transc = 'N';
383 0 : else *(unsigned char *)transc = 'T';
384 :
385 0 : sp_dgemm(transc, "N", A->nrow, nrhs, A->ncol, one, A,
386 : x, ldx, zero, rhs, ldc);
387 :
388 0 : }
389 :
390 : /*! \brief Fills a double precision array with a given value.
391 : */
392 : void
393 0 : dfill(double *a, int alen, double dval)
394 : {
395 : register int i;
396 0 : for (i = 0; i < alen; i++) a[i] = dval;
397 0 : }
398 :
399 :
400 :
401 : /*! \brief Check the inf-norm of the error vector
402 : */
403 0 : 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 0 : Xstore = X->Store;
411 0 : Xmat = Xstore->nzval;
412 :
413 0 : for (j = 0; j < nrhs; j++) {
414 0 : soln_work = &Xmat[j*Xstore->lda];
415 : err = xnorm = 0.0;
416 0 : for (i = 0; i < X->nrow; i++) {
417 0 : err = SUPERLU_MAX(err, fabs(soln_work[i] - xtrue[i]));
418 0 : xnorm = SUPERLU_MAX(xnorm, fabs(soln_work[i]));
419 : }
420 0 : err = err / xnorm;
421 : printf("||X - Xtrue||/||X|| = %e\n", err);
422 : }
423 0 : }
424 :
425 :
426 :
427 : /*! \brief Print performance of the code. */
428 : void
429 0 : dPrintPerf(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage,
430 : double rpg, double rcond, double *ferr,
431 : const double *berr, const char *equed, SuperLUStat_t *stat)
432 : {
433 : SCformat *Lstore;
434 : NCformat *Ustore;
435 : double *utime;
436 : flops_t *ops;
437 :
438 0 : utime = stat->utime;
439 0 : ops = stat->ops;
440 :
441 0 : if ( utime[FACT] != 0. )
442 0 : printf("Factor flops = %e\tMflops = %8.2f\n", ops[FACT],
443 0 : ops[FACT]*1e-6/utime[FACT]);
444 0 : printf("Identify relaxed snodes = %8.2f\n", utime[RELAX]);
445 0 : if ( utime[SOLVE] != 0. )
446 0 : printf("Solve flops = %.0f, Mflops = %8.2f\n", ops[SOLVE],
447 0 : ops[SOLVE]*1e-6/utime[SOLVE]);
448 :
449 0 : Lstore = (SCformat *) L->Store;
450 0 : Ustore = (NCformat *) U->Store;
451 0 : printf("\tNo of nonzeros in factor L = %lld\n", (long long) Lstore->nnz);
452 0 : printf("\tNo of nonzeros in factor U = %lld\n", (long long) Ustore->nnz);
453 0 : printf("\tNo of nonzeros in L+U = %lld\n", (long long) Lstore->nnz + Ustore->nnz);
454 :
455 0 : printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
456 0 : mem_usage->for_lu/1e6, mem_usage->total_needed/1e6);
457 0 : printf("Number of memory expansions: %d\n", stat->expansions);
458 :
459 : printf("\tFactor\tMflops\tSolve\tMflops\tEtree\tEquil\tRcond\tRefine\n");
460 0 : printf("PERF:%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f\n",
461 0 : utime[FACT], ops[FACT]*1e-6/utime[FACT],
462 0 : 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 0 : printf("NUM:\t%e\t%e\t%e\t%e\t%s\n",
467 : rpg, rcond, ferr[0], berr[0], equed);
468 :
469 0 : }
470 :
471 :
472 : int
473 0 : print_double_vec(const char *what, int n, const double *vec)
474 : {
475 : int i;
476 : printf("%s: n %d\n", what, n);
477 0 : for (i = 0; i < n; ++i) printf("%d\t%f\n", i, vec[i]);
478 0 : return 0;
479 : }
|
