1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 #include <stdlib.h>
9 #include <string.h>
10 #include <stdio.h>
11
12 #include "utilities/_hypre_utilities.h"
13 #include "IJ_mv/HYPRE_IJ_mv.h"
14 #include "parcsr_mv/_hypre_parcsr_mv.h"
15 #include "parcsr_ls/_hypre_parcsr_ls.h"
16 #include "seq_mv/seq_mv.h"
17
18 #include "HYPRE_FEI.h"
19
20 extern void hypre_qsort0(int*, int, int);
21 extern void hypre_qsort1(int*, double*, int, int);
22
23 #define habs(x) ((x) > 0.0 ? x : -(x))
24
25 /***************************************************************************/
26 /* reading a matrix from a file in ija format (first row : nrows, nnz) */
27 /* (read by a single processor) */
28 /*-------------------------------------------------------------------------*/
29
HYPRE_LSI_Get_IJAMatrixFromFile(double ** val,int ** ia,int ** ja,int * N,double ** rhs,char * matfile,char * rhsfile)30 void HYPRE_LSI_Get_IJAMatrixFromFile(double **val, int **ia,
31 int **ja, int *N, double **rhs, char *matfile, char *rhsfile)
32 {
33 int i, j, Nrows, nnz, icount, rowindex, colindex, curr_row;
34 int k, m, *mat_ia, *mat_ja, ncnt, rnum;
35 double dtemp, *mat_a, value, *rhs_local;
36 FILE *fp;
37
38 /*------------------------------------------------------------------*/
39 /* read matrix file */
40 /*------------------------------------------------------------------*/
41
42 printf("Reading matrix file = %s \n", matfile );
43 fp = fopen( matfile, "r" );
44 if ( fp == NULL ) {
45 printf("Error : file open error (filename=%s).\n", matfile);
46 exit(1);
47 }
48 fscanf(fp, "%d %d", &Nrows, &nnz);
49 if ( Nrows <= 0 || nnz <= 0 ) {
50 printf("Error : nrows,nnz = %d %d\n", Nrows, nnz);
51 exit(1);
52 }
53 mat_ia = hypre_TAlloc(int, (Nrows+1) , HYPRE_MEMORY_HOST);
54 mat_ja = hypre_TAlloc(int, nnz , HYPRE_MEMORY_HOST);
55 mat_a = hypre_TAlloc(double, nnz , HYPRE_MEMORY_HOST);
56 mat_ia[0] = 0;
57
58 curr_row = 0;
59 icount = 0;
60 for ( i = 0; i < nnz; i++ ) {
61 fscanf(fp, "%d %d %lg", &rowindex, &colindex, &value);
62 rowindex--;
63 colindex--;
64 if ( rowindex != curr_row ) mat_ia[++curr_row] = icount;
65 if ( rowindex < 0 || rowindex >= Nrows )
66 printf("Error reading row %d (curr_row = %d)\n", rowindex, curr_row);
67 if ( colindex < 0 || colindex >= Nrows )
68 printf("Error reading col %d (rowindex = %d)\n", colindex, rowindex);
69 /*if ( value != 0.0 ) {*/
70 mat_ja[icount] = colindex;
71 mat_a[icount++] = value;
72 /*}*/
73 }
74 fclose(fp);
75 for ( i = curr_row+1; i <= Nrows; i++ ) mat_ia[i] = icount;
76 (*val) = mat_a;
77 (*ia) = mat_ia;
78 (*ja) = mat_ja;
79 (*N) = Nrows;
80 printf("matrix has %6d rows and %7d nonzeros\n", Nrows, mat_ia[Nrows]);
81
82 /*------------------------------------------------------------------*/
83 /* read rhs file */
84 /*------------------------------------------------------------------*/
85
86 printf("reading rhs file = %s \n", rhsfile );
87 fp = fopen( rhsfile, "r" );
88 if ( fp == NULL ) {
89 printf("Error : file open error (filename=%s).\n", rhsfile);
90 exit(1);
91 }
92 fscanf(fp, "%d", &ncnt);
93 if ( ncnt <= 0 || ncnt != Nrows) {
94 printf("Error : nrows = %d \n", ncnt);
95 exit(1);
96 }
97 fflush(stdout);
98 rhs_local = hypre_TAlloc(double, Nrows , HYPRE_MEMORY_HOST);
99 m = 0;
100 for ( k = 0; k < ncnt; k++ ) {
101 fscanf(fp, "%d %lg", &rnum, &dtemp);
102 rhs_local[rnum-1] = dtemp; m++;
103 }
104 fflush(stdout);
105 ncnt = m;
106 fclose(fp);
107 (*rhs) = rhs_local;
108 printf("reading rhs done \n");
109 for ( i = 0; i < Nrows; i++ ) {
110 for ( j = mat_ia[i]; j < mat_ia[i+1]; j++ )
111 mat_ja[j]++;
112 }
113 printf("returning from reading matrix\n");
114 }
115
116
117 /***************************************************************************/
118 /* HYPRE_LSI_Search - this is a modification of hypre_BinarySearch */
119 /*-------------------------------------------------------------------------*/
120
HYPRE_LSI_Search(int * list,int value,int list_length)121 int HYPRE_LSI_Search(int *list,int value,int list_length)
122 {
123 int low, high, m;
124 int not_found = 1;
125
126 low = 0;
127 high = list_length-1;
128 while (not_found && low <= high)
129 {
130 m = (low + high) / 2;
131 if (value < list[m])
132 {
133 high = m - 1;
134 }
135 else if (value > list[m])
136 {
137 low = m + 1;
138 }
139 else
140 {
141 not_found = 0;
142 return m;
143 }
144 }
145 return -(low+1);
146 }
147
148 /* ************************************************************************ */
149 /* Given a sorted list of indices and the key, find the position of the */
150 /* key in the list. If not found, return the index of the position */
151 /* corresponding to where it would have been stored. */
152 /* (borrowed from the search routine in ML) */
153 /* ------------------------------------------------------------------------ */
154
HYPRE_LSI_Search2(int key,int nlist,int * list)155 int HYPRE_LSI_Search2(int key, int nlist, int *list)
156 {
157 int nfirst, nlast, nmid, found, index;
158
159 if (nlist <= 0) return -1;
160 nfirst = 0;
161 nlast = nlist-1;
162 if (key > list[nlast]) return -(nlast+1);
163 if (key < list[nfirst]) return -(nfirst+1);
164 found = 0;
165 while ((found == 0) && ((nlast-nfirst)>1)) {
166 nmid = (nfirst + nlast) / 2;
167 if (key == list[nmid]) {index = nmid; found = 1;}
168 else if (key > list[nmid]) nfirst = nmid;
169 else nlast = nmid;
170 }
171 if (found == 1) return index;
172 else if (key == list[nfirst]) return nfirst;
173 else if (key == list[nlast]) return nlast;
174 else return -(nfirst+1);
175 }
176
177 /* ************************************************************************ */
178 /* this function extracts the matrix in a CSR format */
179 /* ------------------------------------------------------------------------ */
180
HYPRE_LSI_GetParCSRMatrix(HYPRE_IJMatrix Amat,int nrows,int nnz,int * ia_ptr,int * ja_ptr,double * a_ptr)181 int HYPRE_LSI_GetParCSRMatrix(HYPRE_IJMatrix Amat, int nrows, int nnz,
182 int *ia_ptr, int *ja_ptr, double *a_ptr)
183 {
184 int nz, i, j, ierr, rowSize, *colInd, nz_ptr, *colInd2;
185 int firstNnz;
186 double *colVal, *colVal2;
187 HYPRE_ParCSRMatrix A_csr;
188
189 nz = 0;
190 nz_ptr = 0;
191 ia_ptr[0] = nz_ptr;
192
193 /* ---old_IJ----------------------------------------------------------- */
194 /*A_csr = (HYPRE_ParCSRMatrix) HYPRE_IJMatrixGetLocalStorage(Amat);*/
195 /* ---new_IJ----------------------------------------------------------- */
196 HYPRE_IJMatrixGetObject(Amat, (void**) &A_csr);
197 /* -------------------------------------------------------------------- */
198
199 for ( i = 0; i < nrows; i++ )
200 {
201 ierr = HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
202 hypre_assert(!ierr);
203 colInd2 = hypre_TAlloc(int, rowSize , HYPRE_MEMORY_HOST);
204 colVal2 = hypre_TAlloc(double, rowSize , HYPRE_MEMORY_HOST);
205 for ( j = 0; j < rowSize; j++ )
206 {
207 colInd2[j] = colInd[j];
208 colVal2[j] = colVal[j];
209 }
210 hypre_qsort1(colInd2, colVal2, 0, rowSize-1);
211 for ( j = 0; j < rowSize-1; j++ )
212 if ( colInd2[j] == colInd2[j+1] )
213 printf("HYPRE_LSI_GetParCSRMatrix-duplicate colind at row %d \n",i);
214
215 firstNnz = 0;
216 for ( j = 0; j < rowSize; j++ )
217 {
218 if ( colVal2[j] != 0.0 )
219 {
220 if (nz_ptr > 0 && firstNnz > 0 && colInd2[j] == ja_ptr[nz_ptr-1])
221 {
222 a_ptr[nz_ptr-1] += colVal2[j];
223 printf("HYPRE_LSI_GetParCSRMatrix:: repeated col in row %d\n",i);
224 }
225 else
226 {
227 ja_ptr[nz_ptr] = colInd2[j];
228 a_ptr[nz_ptr++] = colVal2[j];
229 if ( nz_ptr > nnz )
230 {
231 printf("HYPRE_LSI_GetParCSRMatrix Error (1) - %d %d.\n",i,
232 nrows);
233 exit(1);
234 }
235 firstNnz++;
236 }
237 } else nz++;
238 }
239 hypre_TFree(colInd2, HYPRE_MEMORY_HOST);
240 hypre_TFree(colVal2, HYPRE_MEMORY_HOST);
241 ia_ptr[i+1] = nz_ptr;
242 ierr = HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
243 hypre_assert(!ierr);
244 }
245 /*
246 if ( nnz != nz_ptr )
247 {
248 printf("HYPRE_LSI_GetParCSRMatrix note : matrix sparsity has been \n");
249 printf(" changed since matConfigure - %d > %d ?\n", nnz, nz_ptr);
250 printf(" number of zeros = %d \n", nz );
251 }
252 */
253 return nz_ptr;
254 }
255
256 /* ******************************************************************** */
257 /* sort integers */
258 /* -------------------------------------------------------------------- */
259
HYPRE_LSI_qsort1a(int * ilist,int * ilist2,int left,int right)260 void HYPRE_LSI_qsort1a( int *ilist, int *ilist2, int left, int right)
261 {
262 int i, last, mid, itemp;
263
264 if (left >= right) return;
265 mid = (left + right) / 2;
266 itemp = ilist[left];
267 ilist[left] = ilist[mid];
268 ilist[mid] = itemp;
269 itemp = ilist2[left];
270 ilist2[left] = ilist2[mid];
271 ilist2[mid] = itemp;
272 last = left;
273 for (i = left+1; i <= right; i++)
274 {
275 if (ilist[i] < ilist[left])
276 {
277 last++;
278 itemp = ilist[last];
279 ilist[last] = ilist[i];
280 ilist[i] = itemp;
281 itemp = ilist2[last];
282 ilist2[last] = ilist2[i];
283 ilist2[i] = itemp;
284 }
285 }
286 itemp = ilist[left];
287 ilist[left] = ilist[last];
288 ilist[last] = itemp;
289 itemp = ilist2[left];
290 ilist2[left] = ilist2[last];
291 ilist2[last] = itemp;
292 HYPRE_LSI_qsort1a(ilist, ilist2, left, last-1);
293 HYPRE_LSI_qsort1a(ilist, ilist2, last+1, right);
294 }
295
296 /* ******************************************************************** */
297 /* sort a given list in increasing order */
298 /* -------------------------------------------------------------------- */
299
HYPRE_LSI_SplitDSort2(double * dlist,int nlist,int * ilist,int limit)300 int HYPRE_LSI_SplitDSort2(double *dlist, int nlist, int *ilist, int limit)
301 {
302 int itemp, *iarray1, *iarray2, count1, count2, i;
303 double dtemp, *darray1, *darray2;
304
305 if ( nlist <= 1 ) return 0;
306 if ( nlist == 2 )
307 {
308 if ( dlist[0] < dlist[1] )
309 {
310 dtemp = dlist[0]; dlist[0] = dlist[1]; dlist[1] = dtemp;
311 itemp = ilist[0]; ilist[0] = ilist[1]; ilist[1] = itemp;
312 }
313 return 0;
314 }
315 count1 = 0;
316 count2 = 0;
317 iarray1 = hypre_TAlloc(int, 2 * nlist , HYPRE_MEMORY_HOST);
318 iarray2 = iarray1 + nlist;
319 darray1 = hypre_TAlloc(double, 2 * nlist , HYPRE_MEMORY_HOST);
320 darray2 = darray1 + nlist;
321
322 if ( darray2 == NULL )
323 {
324 printf("ERROR : malloc\n");
325 exit(1);
326 }
327 dtemp = dlist[0];
328 itemp = ilist[0];
329 for ( i = 1; i < nlist; i++ )
330 {
331 if (dlist[i] >= dtemp )
332 {
333 darray1[count1] = dlist[i];
334 iarray1[count1++] = ilist[i];
335 }
336 else
337 {
338 darray2[count2] = dlist[i];
339 iarray2[count2++] = ilist[i];
340 }
341 }
342 dlist[count1] = dtemp;
343 ilist[count1] = itemp;
344 for ( i = 0; i < count1; i++ )
345 {
346 dlist[i] = darray1[i];
347 ilist[i] = iarray1[i];
348 }
349 for ( i = 0; i < count2; i++ )
350 {
351 dlist[count1+1+i] = darray2[i];
352 ilist[count1+1+i] = iarray2[i];
353 }
354 hypre_TFree(darray1, HYPRE_MEMORY_HOST);
355 hypre_TFree(iarray1, HYPRE_MEMORY_HOST);
356 if ( count1+1 == limit ) return 0;
357 else if ( count1+1 < limit )
358 HYPRE_LSI_SplitDSort2(&(dlist[count1+1]),count2,&(ilist[count1+1]),
359 limit-count1-1);
360 else
361 HYPRE_LSI_SplitDSort2( dlist, count1, ilist, limit );
362 return 0;
363 }
364
365 /* ******************************************************************** */
366 /* sort a given list in increasing order */
367 /* -------------------------------------------------------------------- */
368
HYPRE_LSI_SplitDSort(double * dlist,int nlist,int * ilist,int limit)369 int HYPRE_LSI_SplitDSort(double *dlist, int nlist, int *ilist, int limit)
370 {
371 int i, first, last, itemp, cur_index;
372 double dtemp, cur_val;
373
374 if ( nlist <= 1 ) return 0;
375 if ( nlist == 2 )
376 {
377 if ( dlist[0] < dlist[1] )
378 {
379 dtemp = dlist[0]; dlist[0] = dlist[1]; dlist[1] = dtemp;
380 itemp = ilist[0]; ilist[0] = ilist[1]; ilist[1] = itemp;
381 }
382 return 0;
383 }
384
385 first = 0;
386 last = nlist - 1;
387
388 do
389 {
390 cur_index = first;
391 cur_val = dlist[cur_index];
392
393 for ( i = first+1; i <= last; i++ )
394 {
395 if ( dlist[i] > cur_val )
396 {
397 cur_index++;
398 itemp = ilist[cur_index];
399 ilist[cur_index] = ilist[i];
400 ilist[i] = itemp;
401 dtemp = dlist[cur_index];
402 dlist[cur_index] = dlist[i];
403 dlist[i] = dtemp;
404 }
405 }
406 itemp = ilist[cur_index];
407 ilist[cur_index] = ilist[first];
408 ilist[first] = itemp;
409 dtemp = dlist[cur_index];
410 dlist[cur_index] = dlist[first];
411 dlist[first] = dtemp;
412
413 if ( cur_index > limit ) last = cur_index - 1;
414 else if ( cur_index < limit ) first = cur_index + 1;
415 } while ( cur_index != limit );
416
417 return 0;
418 }
419
420 /* ******************************************************************** */
421 /* copy from one vector to another (identity preconditioning) */
422 /* -------------------------------------------------------------------- */
423
HYPRE_LSI_SolveIdentity(HYPRE_Solver solver,HYPRE_ParCSRMatrix Amat,HYPRE_ParVector b,HYPRE_ParVector x)424 int HYPRE_LSI_SolveIdentity(HYPRE_Solver solver, HYPRE_ParCSRMatrix Amat,
425 HYPRE_ParVector b, HYPRE_ParVector x)
426 {
427 (void) solver;
428 (void) Amat;
429 HYPRE_ParVectorCopy( b, x );
430 return 0;
431 }
432
433 /* ******************************************************************** */
434 /* Cuthill McKee reordering algorithm */
435 /* -------------------------------------------------------------------- */
436
HYPRE_LSI_Cuthill(int n,int * ia,int * ja,double * aa,int * order_array,int * reorder_array)437 int HYPRE_LSI_Cuthill(int n, int *ia, int *ja, double *aa, int *order_array,
438 int *reorder_array)
439 {
440 int nnz, *nz_array, cnt, i, j, *tag_array, *queue, nqueue, qhead;
441 int root, norder, mindeg, *ia2, *ja2;
442 double *aa2;
443
444 nz_array = hypre_TAlloc(int, n , HYPRE_MEMORY_HOST);
445 nnz = ia[n];
446 for ( i = 0; i < n; i++ ) nz_array[i] = ia[i+1] - ia[i];
447 tag_array = hypre_TAlloc(int, n , HYPRE_MEMORY_HOST);
448 queue = hypre_TAlloc(int, n , HYPRE_MEMORY_HOST);
449 for ( i = 0; i < n; i++ ) tag_array[i] = 0;
450 norder = 0;
451 mindeg = 10000000;
452 root = -1;
453 for ( i = 0; i < n; i++ )
454 {
455 if ( nz_array[i] == 1 )
456 {
457 tag_array[i] = 1;
458 order_array[norder++] = i;
459 reorder_array[i] = norder-1;
460 }
461 else if ( nz_array[i] < mindeg )
462 {
463 mindeg = nz_array[i];
464 root = i;
465 }
466 }
467 if ( root == -1 )
468 {
469 printf("HYPRE_LSI_Cuthill ERROR : Amat is diagonal\n");
470 exit(1);
471 }
472 nqueue = 0;
473 queue[nqueue++] = root;
474 qhead = 0;
475 tag_array[root] = 1;
476 while ( qhead < nqueue )
477 {
478 root = queue[qhead++];
479 order_array[norder++] = root;
480 reorder_array[root] = norder - 1;
481 for ( j = ia[root]; j < ia[root+1]; j++ )
482 {
483 if ( tag_array[ja[j]] == 0 )
484 {
485 tag_array[ja[j]] = 1;
486 queue[nqueue++] = ja[j];
487 }
488 }
489 if ( qhead == nqueue && norder < n )
490 for ( j = 0; j < n; j++ )
491 if ( tag_array[j] == 0 ) queue[nqueue++] = j;
492 }
493 ia2 = hypre_TAlloc(int, (n+1) , HYPRE_MEMORY_HOST);
494 ja2 = hypre_TAlloc(int, nnz , HYPRE_MEMORY_HOST);
495 aa2 = hypre_TAlloc(double, nnz , HYPRE_MEMORY_HOST);
496 ia2[0] = 0;
497 nnz = 0;
498 for ( i = 0; i < n; i++ )
499 {
500 cnt = order_array[i];
501 for ( j = ia[cnt]; j < ia[cnt+1]; j++ )
502 {
503 ja2[nnz] = ja[j];
504 aa2[nnz++] = aa[j];
505 }
506 ia2[i+1] = nnz;
507 }
508 for ( i = 0; i < nnz; i++ ) ja[i] = reorder_array[ja2[i]];
509 for ( i = 0; i < nnz; i++ ) aa[i] = aa2[i];
510 for ( i = 0; i <= n; i++ ) ia[i] = ia2[i];
511 hypre_TFree(ia2, HYPRE_MEMORY_HOST);
512 hypre_TFree(ja2, HYPRE_MEMORY_HOST);
513 hypre_TFree(aa2, HYPRE_MEMORY_HOST);
514 hypre_TFree(nz_array, HYPRE_MEMORY_HOST);
515 hypre_TFree(tag_array, HYPRE_MEMORY_HOST);
516 hypre_TFree(queue, HYPRE_MEMORY_HOST);
517 return 0;
518 }
519
520 /* ******************************************************************** */
521 /* matrix of a dense matrix */
522 /* -------------------------------------------------------------------- */
523
HYPRE_LSI_MatrixInverse(double ** Amat,int ndim,double *** Cmat)524 int HYPRE_LSI_MatrixInverse( double **Amat, int ndim, double ***Cmat )
525 {
526 int i, j, k;
527 double denom, **Bmat, dmax;
528
529 (*Cmat) = NULL;
530 if ( ndim == 1 )
531 {
532 if ( habs(Amat[0][0]) <= 1.0e-16 ) return -1;
533 Bmat = hypre_TAlloc(double*, ndim , HYPRE_MEMORY_HOST);
534 for ( i = 0; i < ndim; i++ )
535 Bmat[i] = hypre_TAlloc(double, ndim , HYPRE_MEMORY_HOST);
536 Bmat[0][0] = 1.0 / Amat[0][0];
537 (*Cmat) = Bmat;
538 return 0;
539 }
540 if ( ndim == 2 )
541 {
542 denom = Amat[0][0] * Amat[1][1] - Amat[0][1] * Amat[1][0];
543 if ( habs( denom ) <= 1.0e-16 ) return -1;
544 Bmat = hypre_TAlloc(double*, ndim , HYPRE_MEMORY_HOST);
545 for ( i = 0; i < ndim; i++ )
546 Bmat[i] = hypre_TAlloc(double, ndim , HYPRE_MEMORY_HOST);
547 Bmat[0][0] = Amat[1][1] / denom;
548 Bmat[1][1] = Amat[0][0] / denom;
549 Bmat[0][1] = - ( Amat[0][1] / denom );
550 Bmat[1][0] = - ( Amat[1][0] / denom );
551 (*Cmat) = Bmat;
552 return 0;
553 }
554 else
555 {
556 Bmat = hypre_TAlloc(double*, ndim , HYPRE_MEMORY_HOST);
557 for ( i = 0; i < ndim; i++ )
558 {
559 Bmat[i] = hypre_TAlloc(double, ndim , HYPRE_MEMORY_HOST);
560 for ( j = 0; j < ndim; j++ ) Bmat[i][j] = 0.0;
561 Bmat[i][i] = 1.0;
562 }
563 for ( i = 1; i < ndim; i++ )
564 {
565 for ( j = 0; j < i; j++ )
566 {
567 if ( habs(Amat[j][j]) < 1.0e-16 ) return -1;
568 denom = Amat[i][j] / Amat[j][j];
569 for ( k = 0; k < ndim; k++ )
570 {
571 Amat[i][k] -= denom * Amat[j][k];
572 Bmat[i][k] -= denom * Bmat[j][k];
573 }
574 }
575 }
576 for ( i = ndim-2; i >= 0; i-- )
577 {
578 for ( j = ndim-1; j >= i+1; j-- )
579 {
580 if ( habs(Amat[j][j]) < 1.0e-16 ) return -1;
581 denom = Amat[i][j] / Amat[j][j];
582 for ( k = 0; k < ndim; k++ )
583 {
584 Amat[i][k] -= denom * Amat[j][k];
585 Bmat[i][k] -= denom * Bmat[j][k];
586 }
587 }
588 }
589 for ( i = 0; i < ndim; i++ )
590 {
591 denom = Amat[i][i];
592 if ( habs(denom) < 1.0e-16 ) return -1;
593 for ( j = 0; j < ndim; j++ ) Bmat[i][j] /= denom;
594 }
595
596 for ( i = 0; i < ndim; i++ )
597 for ( j = 0; j < ndim; j++ )
598 if ( habs(Bmat[i][j]) < 1.0e-17 ) Bmat[i][j] = 0.0;
599 dmax = 0.0;
600 for ( i = 0; i < ndim; i++ )
601 {
602 for ( j = 0; j < ndim; j++ )
603 if ( habs(Bmat[i][j]) > dmax ) dmax = habs(Bmat[i][j]);
604 /*
605 for ( j = 0; j < ndim; j++ )
606 if ( habs(Bmat[i][j]/dmax) < 1.0e-15 ) Bmat[i][j] = 0.0;
607 */
608 }
609 (*Cmat) = Bmat;
610 if ( dmax > 1.0e6 ) return 1;
611 else return 0;
612 }
613 }
614
615 /* ******************************************************************** */
616 /* find the separators of a matrix */
617 /* -------------------------------------------------------------------- */
618
HYPRE_LSI_PartitionMatrix(int nRows,int startRow,int * rowLengths,int ** colIndices,double ** colValues,int * nLabels,int ** labels)619 int HYPRE_LSI_PartitionMatrix( int nRows, int startRow, int *rowLengths,
620 int **colIndices, double **colValues,
621 int *nLabels, int **labels)
622 {
623 int irow, rowCnt, labelNum, *localLabels, actualNRows;
624 int jcol, root, indHead, indTail, *indSet, index;
625
626 /*----------------------------------------------------------------*/
627 /* search for constraint rows */
628 /*----------------------------------------------------------------*/
629
630 for ( irow = nRows-1; irow >= 0; irow-- )
631 {
632 index = irow + startRow;
633 for ( jcol = 0; jcol < rowLengths[irow]; jcol++ )
634 if (colIndices[irow][jcol] == index && colValues[irow][jcol] != 0.0)
635 break;
636 if ( jcol != rowLengths[irow] ) break;
637 }
638 (*nLabels) = actualNRows = irow + 1;
639
640 /*----------------------------------------------------------------*/
641 /* search for constraint rows */
642 /*----------------------------------------------------------------*/
643
644 localLabels = hypre_TAlloc(int, actualNRows , HYPRE_MEMORY_HOST);
645 for ( irow = 0; irow < actualNRows; irow++ ) localLabels[irow] = -1;
646 indSet = hypre_TAlloc(int, actualNRows , HYPRE_MEMORY_HOST);
647
648 labelNum = 0;
649 rowCnt = actualNRows;
650
651 while ( rowCnt > 0 )
652 {
653 root = -1;
654 for ( irow = 0; irow < actualNRows; irow++ )
655 if ( localLabels[irow] == -1 ) {root = irow; break;}
656 if ( root == -1 )
657 {
658 printf("HYPRE_LSI_PartitionMatrix : something wrong.\n");
659 exit(1);
660 }
661 indHead = indTail = 0;
662 localLabels[root] = labelNum;
663 rowCnt--;
664 for ( jcol = 0; jcol < rowLengths[root]; jcol++ )
665 {
666 index = colIndices[root][jcol] - startRow;
667 if ( index >= 0 && index < actualNRows && localLabels[index] < 0 )
668 {
669 indSet[indTail++] = index;
670 localLabels[index] = labelNum;
671 }
672 }
673 while ( (indTail - indHead) > 0 )
674 {
675 root = indSet[indHead++];
676 rowCnt--;
677 for ( jcol = 0; jcol < rowLengths[root]; jcol++ )
678 {
679 index = colIndices[root][jcol] - startRow;
680 if ( index >= 0 && index < actualNRows && localLabels[index] < 0 )
681 {
682 indSet[indTail++] = index;
683 localLabels[index] = labelNum;
684 }
685 }
686 }
687 labelNum++;
688 }
689 if ( labelNum > 4 )
690 {
691 printf("HYPRE_LSI_PartitionMatrix : number of labels %d too large.\n",
692 labelNum+1);
693 hypre_TFree(localLabels, HYPRE_MEMORY_HOST);
694 (*nLabels) = 0;
695 (*labels) = NULL;
696 }
697 else
698 {
699 printf("HYPRE_LSI_PartitionMatrix : number of labels = %d.\n",
700 labelNum);
701 (*labels) = localLabels;
702 }
703 hypre_TFree(indSet, HYPRE_MEMORY_HOST);
704 return 0;
705 }
706
707