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 /******************************************************************************
9 *
10 * HYPRE_DDICT interface
11 *
12 *****************************************************************************/
13
14 #include <stdlib.h>
15 #include <stdio.h>
16 #include <math.h>
17
18 #include "utilities/_hypre_utilities.h"
19 #include "HYPRE.h"
20 #include "IJ_mv/HYPRE_IJ_mv.h"
21 #include "parcsr_mv/HYPRE_parcsr_mv.h"
22 #include "parcsr_mv/_hypre_parcsr_mv.h"
23 #include "parcsr_ls/HYPRE_parcsr_ls.h"
24 #include "HYPRE_MHMatrix.h"
25
26 #ifdef HAVE_ML
27
28 #include "ml_struct.h"
29 #include "ml_aggregate.h"
30
31 #endif
32
33 #include "HYPRE_MHMatrix.h"
34 #include "HYPRE_FEI.h"
35 typedef struct HYPRE_LSI_DDICT_Struct
36 {
37 MPI_Comm comm;
38 MH_Matrix *mh_mat;
39 double thresh;
40 double fillin;
41 int Nrows;
42 int extNrows;
43 int *mat_ja;
44 double *mat_aa;
45 int outputLevel;
46 }
47 HYPRE_LSI_DDICT;
48
49 extern int HYPRE_LSI_MLConstructMHMatrix(HYPRE_ParCSRMatrix,MH_Matrix *,
50 MPI_Comm, int *, MH_Context *);
51 extern int HYPRE_LSI_DDICTComposeOverlappedMatrix(MH_Matrix *, int *,
52 int **recv_lengths, int **int_buf, double **dble_buf,
53 int **sindex_array, int **sindex_array2, int *offset);
54 extern int HYPRE_LSI_DDICTGetRowLengths(MH_Matrix *Amat, int *leng, int **);
55 extern int HYPRE_LIS_DDICTGetOffProcRows(MH_Matrix *Amat, int leng, int *,
56 int Noffset, int *map, int *map2, int **int_buf,
57 double **dble_buf);
58 extern int HYPRE_LSI_DDICTDecompose(HYPRE_LSI_DDICT *ict_ptr,MH_Matrix *Amat,
59 int total_recv_leng, int *recv_lengths, int *ext_ja,
60 double *ext_aa, int *map, int *map2, int Noffset);
61 extern void HYPRE_LSI_qsort1a(int *, int *, int, int);
62 extern int HYPRE_LSI_SplitDSort(double *,int,int*,int);
63 extern int HYPRE_LSI_Search(int *, int, int);
64
65 extern int HYPRE_LSI_DDICTFactorize(HYPRE_LSI_DDICT *ict_ptr, double *mat_aa,
66 int *mat_ja, int *mat_ia, double *rowNorms);
67
68 extern int MH_ExchBdry(double *, void *);
69 extern int MH_ExchBdryBack(double *, void *, int *, double **, int **);
70 extern int MH_GetRow(void *, int, int *, int, int *, double *, int *);
71
72 #define habs(x) ((x) > 0 ? (x) : -(x))
73
74 /*****************************************************************************/
75 /* HYPRE_LSI_DDICTCreate - Return a DDICT preconditioner object "solver". */
76 /*---------------------------------------------------------------------------*/
77
HYPRE_LSI_DDICTCreate(MPI_Comm comm,HYPRE_Solver * solver)78 int HYPRE_LSI_DDICTCreate( MPI_Comm comm, HYPRE_Solver *solver )
79 {
80 HYPRE_LSI_DDICT *ict_ptr;
81
82 ict_ptr = hypre_TAlloc(HYPRE_LSI_DDICT, 1, HYPRE_MEMORY_HOST);
83
84 if (ict_ptr == NULL) return 1;
85
86 ict_ptr->comm = comm;
87 ict_ptr->mh_mat = NULL;
88 ict_ptr->fillin = 0.0;
89 ict_ptr->thresh = 0.0; /* defaults */
90 ict_ptr->mat_ja = NULL;
91 ict_ptr->mat_aa = NULL;
92 ict_ptr->outputLevel = 0;
93
94 *solver = (HYPRE_Solver) ict_ptr;
95
96 return 0;
97 }
98
99 /*****************************************************************************/
100 /* HYPRE_LSI_DDICTDestroy - Destroy a DDICT object. */
101 /*---------------------------------------------------------------------------*/
102
HYPRE_LSI_DDICTDestroy(HYPRE_Solver solver)103 int HYPRE_LSI_DDICTDestroy( HYPRE_Solver solver )
104 {
105 int i;
106 HYPRE_LSI_DDICT *ict_ptr;
107
108 ict_ptr = (HYPRE_LSI_DDICT *) solver;
109 hypre_TFree(ict_ptr->mat_ja, HYPRE_MEMORY_HOST);
110 hypre_TFree(ict_ptr->mat_aa, HYPRE_MEMORY_HOST);
111 if ( ict_ptr->mh_mat != NULL )
112 {
113 hypre_TFree(ict_ptr->mh_mat->sendProc, HYPRE_MEMORY_HOST);
114 hypre_TFree(ict_ptr->mh_mat->sendLeng, HYPRE_MEMORY_HOST);
115 hypre_TFree(ict_ptr->mh_mat->recvProc, HYPRE_MEMORY_HOST);
116 hypre_TFree(ict_ptr->mh_mat->recvLeng, HYPRE_MEMORY_HOST);
117 for ( i = 0; i < ict_ptr->mh_mat->sendProcCnt; i++ )
118 hypre_TFree(ict_ptr->mh_mat->sendList[i], HYPRE_MEMORY_HOST);
119 hypre_TFree(ict_ptr->mh_mat->sendList, HYPRE_MEMORY_HOST);
120 hypre_TFree(ict_ptr, HYPRE_MEMORY_HOST);
121 }
122 ict_ptr->mh_mat = NULL;
123 hypre_TFree(ict_ptr, HYPRE_MEMORY_HOST);
124
125 return 0;
126 }
127
128 /*****************************************************************************/
129 /* HYPRE_LSI_DDICTSetFillin - Set the fill-in parameter. */
130 /*---------------------------------------------------------------------------*/
131
HYPRE_LSI_DDICTSetFillin(HYPRE_Solver solver,double fillin)132 int HYPRE_LSI_DDICTSetFillin(HYPRE_Solver solver, double fillin)
133 {
134 HYPRE_LSI_DDICT *ict_ptr = (HYPRE_LSI_DDICT *) solver;
135
136 ict_ptr->fillin = fillin;
137
138 return 0;
139 }
140
141 /*****************************************************************************/
142 /* HYPRE_LSI_DDICTSetDropTolerance - Set the threshold for dropping */
143 /*---------------------------------------------------------------------------*/
144
HYPRE_LSI_DDICTSetDropTolerance(HYPRE_Solver solver,double thresh)145 int HYPRE_LSI_DDICTSetDropTolerance(HYPRE_Solver solver, double thresh)
146 {
147 HYPRE_LSI_DDICT *ict_ptr = (HYPRE_LSI_DDICT *) solver;
148
149 ict_ptr->thresh = thresh;
150
151 return 0;
152 }
153
154 /*****************************************************************************/
155 /* HYPRE_LSI_DDICTSetOutputLevel - Set debug level */
156 /*---------------------------------------------------------------------------*/
157
HYPRE_LSI_DDICTSetOutputLevel(HYPRE_Solver solver,int level)158 int HYPRE_LSI_DDICTSetOutputLevel(HYPRE_Solver solver, int level)
159 {
160 HYPRE_LSI_DDICT *ict_ptr = (HYPRE_LSI_DDICT *) solver;
161
162 ict_ptr->outputLevel = level;
163
164 return 0;
165 }
166
167 /*****************************************************************************/
168 /* HYPRE_LSI_DDICTSolve - Solve function for DDICT. */
169 /*---------------------------------------------------------------------------*/
170
HYPRE_LSI_DDICTSolve(HYPRE_Solver solver,HYPRE_ParCSRMatrix A,HYPRE_ParVector b,HYPRE_ParVector x)171 int HYPRE_LSI_DDICTSolve( HYPRE_Solver solver, HYPRE_ParCSRMatrix A,
172 HYPRE_ParVector b, HYPRE_ParVector x )
173 {
174 int i, j, Nrows, extNrows, *mat_ja, *ibuf, length;
175 double *rhs, *soln, *dbuf, *mat_aa, *dbuf2, dtmp;
176 HYPRE_LSI_DDICT *ict_ptr = (HYPRE_LSI_DDICT *) solver;
177 MH_Context *context;
178
179 rhs = hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *) b));
180 soln = hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *) x));
181
182 Nrows = ict_ptr->Nrows;
183 extNrows = ict_ptr->extNrows;
184 mat_ja = ict_ptr->mat_ja;
185 mat_aa = ict_ptr->mat_aa;
186
187 if ( extNrows > 0 )
188 {
189 dbuf = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
190 dbuf2 = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
191 for ( i = 0; i < Nrows; i++ ) dbuf[i] = rhs[i];
192 }
193 else dbuf = dbuf2 = NULL;
194
195 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
196 context->Amat = ict_ptr->mh_mat;
197 context->comm = MPI_COMM_WORLD;
198
199 MH_ExchBdry(dbuf, context);
200
201 for ( i = 0; i < extNrows; i++ )
202 {
203 dtmp = dbuf[i];
204 for ( j = mat_ja[i]; j < mat_ja[i+1]; j++ )
205 dtmp -= ( mat_aa[j] * dbuf2[mat_ja[j]] );
206 dbuf2[i] = dtmp * mat_aa[i];
207 }
208 for ( i = extNrows-1; i >= 0; i-- )
209 {
210 dbuf2[i] *= mat_aa[i];
211 dtmp = dbuf2[i];
212 for ( j = mat_ja[i]; j < mat_ja[i+1]; j++ )
213 dbuf2[mat_ja[j]] -= ( dtmp * mat_aa[j] );
214 }
215 hypre_TFree(dbuf, HYPRE_MEMORY_HOST);
216
217 for ( i = 0; i < Nrows; i++ ) soln[i] = dbuf2[i];
218
219 MH_ExchBdryBack(dbuf2, context, &length, &dbuf, &ibuf);
220
221 for ( i = 0; i < length; i++ ) soln[ibuf[i]] = soln[ibuf[i]] + dbuf[i];
222
223 hypre_TFree(ibuf, HYPRE_MEMORY_HOST);
224 hypre_TFree(dbuf, HYPRE_MEMORY_HOST);
225 hypre_TFree(dbuf2, HYPRE_MEMORY_HOST);
226 hypre_TFree(context, HYPRE_MEMORY_HOST);
227
228 return 0;
229 }
230
231 /*****************************************************************************/
232 /* HYPRE_LSI_DDICTSetup - Set up function for LSI_DDICT. */
233 /*---------------------------------------------------------------------------*/
234
HYPRE_LSI_DDICTSetup(HYPRE_Solver solver,HYPRE_ParCSRMatrix A_csr,HYPRE_ParVector b,HYPRE_ParVector x)235 int HYPRE_LSI_DDICTSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix A_csr,
236 HYPRE_ParVector b, HYPRE_ParVector x )
237 {
238 int i, j, offset, total_recv_leng, *recv_lengths=NULL;
239 int *int_buf=NULL, mypid, nprocs, overlap_flag=1,*parray;
240 int *map=NULL, *map2=NULL, *row_partition=NULL,*parray2;
241 double *dble_buf=NULL;
242 HYPRE_LSI_DDICT *ict_ptr = (HYPRE_LSI_DDICT *) solver;
243 MH_Context *context=NULL;
244 MH_Matrix *mh_mat=NULL;
245
246 /* ---------------------------------------------------------------- */
247 /* get the row information in my processors */
248 /* ---------------------------------------------------------------- */
249
250 MPI_Comm_rank(MPI_COMM_WORLD, &mypid);
251 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
252 HYPRE_ParCSRMatrixGetRowPartitioning(A_csr, &row_partition);
253
254 /* ---------------------------------------------------------------- */
255 /* convert the incoming CSR matrix into a MH matrix */
256 /* ---------------------------------------------------------------- */
257
258 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
259 context->comm = MPI_COMM_WORLD;
260 context->globalEqns = row_partition[nprocs];
261 context->partition = hypre_TAlloc(int, (nprocs+1), HYPRE_MEMORY_HOST);
262 for (i=0; i<=nprocs; i++) context->partition[i] = row_partition[i];
263 hypre_TFree( row_partition , HYPRE_MEMORY_HOST);
264 mh_mat = hypre_TAlloc( MH_Matrix, 1, HYPRE_MEMORY_HOST);
265 context->Amat = mh_mat;
266 HYPRE_LSI_MLConstructMHMatrix(A_csr,mh_mat,MPI_COMM_WORLD,
267 context->partition,context);
268
269 /* ---------------------------------------------------------------- */
270 /* compose the enlarged overlapped local matrix */
271 /* ---------------------------------------------------------------- */
272
273 if ( overlap_flag )
274 {
275 HYPRE_LSI_DDICTComposeOverlappedMatrix(mh_mat, &total_recv_leng,
276 &recv_lengths, &int_buf, &dble_buf, &map, &map2,&offset);
277 }
278 else
279 {
280 total_recv_leng = 0;
281 recv_lengths = NULL;
282 int_buf = NULL;
283 dble_buf = NULL;
284 map = NULL;
285 map2 = NULL;
286 parray = hypre_TAlloc(int, nprocs , HYPRE_MEMORY_HOST);
287 parray2 = hypre_TAlloc(int, nprocs , HYPRE_MEMORY_HOST);
288 for ( i = 0; i < nprocs; i++ ) parray2[i] = 0;
289 parray2[mypid] = mh_mat->Nrows;
290 MPI_Allreduce(parray2,parray,nprocs,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
291 offset = 0;
292 for (i = 0; i < mypid; i++) offset += parray[i];
293 hypre_TFree(parray, HYPRE_MEMORY_HOST);
294 hypre_TFree(parray2, HYPRE_MEMORY_HOST);
295 }
296
297 /* ---------------------------------------------------------------- */
298 /* perform ICT decomposition on local matrix */
299 /* ---------------------------------------------------------------- */
300
301 HYPRE_LSI_DDICTDecompose(ict_ptr,mh_mat,total_recv_leng,recv_lengths,
302 int_buf, dble_buf, map,map2, offset);
303
304 if ( mypid == 0 && ict_ptr->outputLevel > 2 )
305 {
306 for ( i = 0; i < ict_ptr->extNrows; i++ )
307 for ( j = ict_ptr->mat_ja[i]; j < ict_ptr->mat_ja[i+1]; j++ )
308 printf("LA(%d,%d) = %e;\n", i+1, ict_ptr->mat_ja[j]+1,
309 ict_ptr->mat_aa[j]);
310 }
311 ict_ptr->mh_mat = mh_mat;
312 hypre_TFree(recv_lengths, HYPRE_MEMORY_HOST);
313 hypre_TFree(int_buf, HYPRE_MEMORY_HOST);
314 hypre_TFree(dble_buf, HYPRE_MEMORY_HOST);
315 hypre_TFree(map, HYPRE_MEMORY_HOST);
316 hypre_TFree(map2, HYPRE_MEMORY_HOST);
317 hypre_TFree(context->partition, HYPRE_MEMORY_HOST);
318 hypre_TFree(context, HYPRE_MEMORY_HOST);
319 return 0;
320 }
321
322 /*****************************************************************************/
323 /* subroutines used for constructing overlapped matrix */
324 /*---------------------------------------------------------------------------*/
325
HYPRE_LSI_DDICTGetRowLengths(MH_Matrix * Amat,int * leng,int ** recv_leng)326 int HYPRE_LSI_DDICTGetRowLengths(MH_Matrix *Amat, int *leng, int **recv_leng)
327 {
328 int i, j, m, mypid, index, *temp_list, allocated_space, length;
329 int nRecv, *recvProc, *recvLeng, *cols, total_recv, mtype, msgtype;
330 int nSend, *sendProc, *sendLeng, **sendList, proc_id, offset;
331 double *vals;
332 MPI_Request *Request;
333 MPI_Status status;
334 MH_Context *context;
335
336 /* ---------------------------------------------------------------- */
337 /* fetch communication information */
338 /* ---------------------------------------------------------------- */
339
340 MPI_Comm_rank(MPI_COMM_WORLD, &mypid);
341 nRecv = Amat->recvProcCnt;
342 nSend = Amat->sendProcCnt;
343 recvProc = Amat->recvProc;
344 recvLeng = Amat->recvLeng;
345 sendProc = Amat->sendProc;
346 sendLeng = Amat->sendLeng;
347 sendList = Amat->sendList;
348 total_recv = 0;
349 for ( i = 0; i < nRecv; i++ ) total_recv += recvLeng[i];
350
351 (*leng) = total_recv;
352 if ( nRecv <= 0 ) (*recv_leng) = NULL;
353
354 MPI_Barrier(MPI_COMM_WORLD);
355
356 /* ---------------------------------------------------------------- */
357 /* post receives for all messages */
358 /* ---------------------------------------------------------------- */
359
360 (*recv_leng) = hypre_TAlloc(int, total_recv , HYPRE_MEMORY_HOST);
361 if (nRecv > 0) Request = hypre_TAlloc(MPI_Request, nRecv, HYPRE_MEMORY_HOST);
362 offset = 0;
363 mtype = 2001;
364 for (i = 0; i < nRecv; i++)
365 {
366 proc_id = recvProc[i];
367 msgtype = mtype;
368 length = recvLeng[i];
369 MPI_Irecv((void *) &((*recv_leng)[offset]), length, MPI_INT, proc_id,
370 msgtype, MPI_COMM_WORLD, &Request[i]);
371 offset += length;
372 }
373
374 /* ---------------------------------------------------------------- */
375 /* write out all messages */
376 /* ---------------------------------------------------------------- */
377
378 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
379 context->Amat = Amat;
380 allocated_space = 100;
381 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
382 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
383 for (i = 0; i < nSend; i++)
384 {
385 proc_id = sendProc[i];
386 length = sendLeng[i];
387 temp_list = hypre_TAlloc(int, sendLeng[i] , HYPRE_MEMORY_HOST);
388 for (j = 0; j < length; j++)
389 {
390 index = sendList[i][j];
391 while (MH_GetRow(context,1,&index,allocated_space,cols,vals,&m)==0)
392 {
393 hypre_TFree(cols, HYPRE_MEMORY_HOST);
394 hypre_TFree(vals, HYPRE_MEMORY_HOST);
395 allocated_space += 200 + 1;
396 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
397 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
398 }
399 temp_list[j] = m;
400 }
401 msgtype = mtype;
402 MPI_Send((void*)temp_list,length,MPI_INT,proc_id,msgtype,MPI_COMM_WORLD);
403 hypre_TFree(temp_list, HYPRE_MEMORY_HOST);
404 }
405 hypre_TFree(cols, HYPRE_MEMORY_HOST);
406 hypre_TFree(vals, HYPRE_MEMORY_HOST);
407 hypre_TFree(context, HYPRE_MEMORY_HOST);
408
409 /* ---------------------------------------------------------------- */
410 /* wait for messages */
411 /* ---------------------------------------------------------------- */
412
413 for ( i = 0; i < nRecv; i++ )
414 {
415 MPI_Wait( &Request[i], &status );
416 }
417
418 if (nRecv > 0)
419 hypre_TFree(Request, HYPRE_MEMORY_HOST);
420 return 0;
421 }
422
423 /*****************************************************************************/
424 /* needed for overlapped smoothers */
425 /*---------------------------------------------------------------------------*/
426
HYPRE_LSI_DDICTGetOffProcRows(MH_Matrix * Amat,int leng,int * recv_leng,int Noffset,int * map,int * map2,int ** int_buf,double ** dble_buf)427 int HYPRE_LSI_DDICTGetOffProcRows(MH_Matrix *Amat, int leng, int *recv_leng,
428 int Noffset, int *map, int *map2, int **int_buf,
429 double **dble_buf)
430 {
431 int i, j, k, m, length, offset, allocated_space, proc_id;
432 int nRecv, nSend, *recvProc, *sendProc, total_recv, mtype, msgtype;
433 int *sendLeng, *recvLeng, **sendList, *cols, *isend_buf, Nrows;
434 int nnz, nnz_offset, index, mypid;
435 double *vals, *send_buf;
436 MPI_Request *request;
437 MPI_Status status;
438 MH_Context *context;
439
440 /* ---------------------------------------------------------------- */
441 /* fetch communication information */
442 /* ---------------------------------------------------------------- */
443
444 MPI_Comm_rank(MPI_COMM_WORLD, &mypid);
445 Nrows = Amat->Nrows;
446 nRecv = Amat->recvProcCnt;
447 nSend = Amat->sendProcCnt;
448 recvProc = Amat->recvProc;
449 recvLeng = Amat->recvLeng;
450 sendProc = Amat->sendProc;
451 sendLeng = Amat->sendLeng;
452 sendList = Amat->sendList;
453 if ( nRecv <= 0 ) { (*int_buf) = NULL; (*dble_buf) = NULL;}
454 total_recv = 0;
455 for ( i = 0; i < leng; i++ ) total_recv += recv_leng[i];
456
457 /* ---------------------------------------------------------------- */
458 /* allocate buffer space */
459 /* ---------------------------------------------------------------- */
460
461 if ( nRecv > 0 )
462 request = hypre_TAlloc(MPI_Request , nRecv, HYPRE_MEMORY_HOST);
463 else request = NULL;
464
465 if ( total_recv > 0 )
466 {
467 (*int_buf) = hypre_TAlloc(int, total_recv , HYPRE_MEMORY_HOST);
468 (*dble_buf) = hypre_TAlloc(double, total_recv , HYPRE_MEMORY_HOST);
469 }
470
471 /* ---------------------------------------------------------------- */
472 /* post receives for all messages */
473 /* ---------------------------------------------------------------- */
474
475 offset = 0;
476 mtype = 2002;
477 nnz_offset = 0;
478 for (i = 0; i < nRecv; i++)
479 {
480 proc_id = recvProc[i];
481 msgtype = mtype;
482 length = recvLeng[i];
483 nnz = 0;
484 for (j = 0; j < length; j++) nnz += recv_leng[offset+j];
485
486 MPI_Irecv((void *) &((*dble_buf)[nnz_offset]), nnz, MPI_DOUBLE,
487 proc_id, msgtype, MPI_COMM_WORLD, request+i);
488 offset += length;
489 nnz_offset += nnz;
490 }
491
492 /* ---------------------------------------------------------------- */
493 /* send rows to other processors */
494 /* ---------------------------------------------------------------- */
495
496 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
497 context->Amat = Amat;
498 mtype = 2002;
499 allocated_space = 100;
500 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
501 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
502 for (i = 0; i < nSend; i++)
503 {
504 proc_id = sendProc[i];
505 length = sendLeng[i];
506 nnz = 0;
507 for (j = 0; j < length; j++)
508 {
509 index = sendList[i][j];
510 while (MH_GetRow(context,1,&index,allocated_space,cols,vals,&m)==0)
511 {
512 hypre_TFree(cols, HYPRE_MEMORY_HOST);
513 hypre_TFree(vals, HYPRE_MEMORY_HOST);
514 allocated_space += 200 + 1;
515 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
516 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
517 }
518 nnz += m;
519 }
520 if ( nnz > 0 ) send_buf = hypre_TAlloc(double, nnz , HYPRE_MEMORY_HOST);
521 offset = 0;
522 for (j = 0; j < length; j++)
523 {
524 index = sendList[i][j];
525 MH_GetRow(context,1,&index,allocated_space,cols,vals,&m);
526 for (k = 0; k < m; k++) send_buf[offset+k] = vals[k];
527 offset += m;
528 }
529 msgtype = mtype;
530 MPI_Send((void*) send_buf, nnz, MPI_DOUBLE, proc_id, msgtype,
531 MPI_COMM_WORLD);
532 if ( nnz > 0 )
533 hypre_TFree(send_buf, HYPRE_MEMORY_HOST);
534 }
535 hypre_TFree(cols, HYPRE_MEMORY_HOST);
536 hypre_TFree(vals, HYPRE_MEMORY_HOST);
537
538 /* ---------------------------------------------------------------- */
539 /* wait for all messages */
540 /* ---------------------------------------------------------------- */
541
542 for (i = 0; i < nRecv; i++)
543 {
544 MPI_Wait(request+i, &status);
545 }
546
547 /* ----------------------------------------------------------- */
548 /* post receives for all messages */
549 /* ----------------------------------------------------------- */
550
551 mtype = 2003;
552 offset = 0;
553 nnz_offset = 0;
554 for (i = 0; i < nRecv; i++)
555 {
556 proc_id = recvProc[i];
557 msgtype = mtype;
558 length = recvLeng[i];
559 nnz = 0;
560 for (j = 0; j < length; j++) nnz += recv_leng[offset+j];
561 MPI_Irecv((void *) &((*int_buf)[nnz_offset]), nnz, MPI_INT,
562 proc_id, msgtype, MPI_COMM_WORLD, request+i);
563 offset += length;
564 nnz_offset += nnz;
565 }
566
567 /* ---------------------------------------------------------------- */
568 /* send rows to other processors */
569 /* ---------------------------------------------------------------- */
570
571 mtype = 2003;
572 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
573 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
574 for (i = 0; i < nSend; i++)
575 {
576 proc_id = sendProc[i];
577 length = sendLeng[i];
578 nnz = 0;
579 for (j = 0; j < length; j++)
580 {
581 index = sendList[i][j];
582 MH_GetRow(context,1,&index,allocated_space,cols,vals,&m);
583 nnz += m;
584 }
585 if ( nnz > 0 ) isend_buf = hypre_TAlloc(int, nnz , HYPRE_MEMORY_HOST);
586 offset = 0;
587 for (j = 0; j < length; j++)
588 {
589 index = sendList[i][j];
590 MH_GetRow(context,1,&index,allocated_space,cols,vals,&m);
591 for (k = 0; k < m; k++)
592 {
593 if ( cols[k] < Nrows ) isend_buf[offset+k] = cols[k] + Noffset;
594 else isend_buf[offset+k] = map[cols[k]-Nrows];
595 }
596 offset += m;
597 }
598 msgtype = mtype;
599 MPI_Send((void*) isend_buf, nnz, MPI_INT, proc_id, msgtype,
600 MPI_COMM_WORLD);
601 if ( nnz > 0 )
602 hypre_TFree(isend_buf, HYPRE_MEMORY_HOST);
603 }
604 hypre_TFree(cols, HYPRE_MEMORY_HOST);
605 hypre_TFree(vals, HYPRE_MEMORY_HOST);
606
607 /* ----------------------------------------------------------- */
608 /* post receives for all messages */
609 /* ----------------------------------------------------------- */
610
611 for (i = 0; i < nRecv; i++)
612 {
613 MPI_Wait(request+i, &status);
614 }
615
616 hypre_TFree(request, HYPRE_MEMORY_HOST);
617 hypre_TFree(context, HYPRE_MEMORY_HOST);
618 return 0;
619 }
620
621 /*****************************************************************************/
622 /* construct an enlarged overlapped local matrix */
623 /*---------------------------------------------------------------------------*/
624
HYPRE_LSI_DDICTComposeOverlappedMatrix(MH_Matrix * mh_mat,int * total_recv_leng,int ** recv_lengths,int ** int_buf,double ** dble_buf,int ** sindex_array,int ** sindex_array2,int * offset)625 int HYPRE_LSI_DDICTComposeOverlappedMatrix(MH_Matrix *mh_mat,
626 int *total_recv_leng, int **recv_lengths, int **int_buf,
627 double **dble_buf, int **sindex_array, int **sindex_array2,
628 int *offset)
629 {
630 int i, nprocs, mypid, Nrows, *proc_array, *proc_array2;
631 int extNrows, NrowsOffset, *index_array, *index_array2;
632 int nRecv, *recvLeng;
633 double *dble_array;
634 MH_Context *context;
635
636 /* ---------------------------------------------------------------- */
637 /* fetch communication information */
638 /* ---------------------------------------------------------------- */
639
640 MPI_Comm_rank(MPI_COMM_WORLD, &mypid);
641 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
642
643 /* ---------------------------------------------------------------- */
644 /* fetch matrix information */
645 /* ---------------------------------------------------------------- */
646
647 nRecv = mh_mat->recvProcCnt;
648 recvLeng = mh_mat->recvLeng;
649 Nrows = mh_mat->Nrows;
650
651 /* ---------------------------------------------------------------- */
652 /* compute the enlarged matrix size */
653 /* ---------------------------------------------------------------- */
654
655 (*total_recv_leng) = 0;
656 for ( i = 0; i < nRecv; i++ ) (*total_recv_leng) += recvLeng[i];
657 extNrows = Nrows + (*total_recv_leng);
658
659 /* ---------------------------------------------------------------- */
660 /* compose NrowsOffset and processor offsets proc_array */
661 /* ---------------------------------------------------------------- */
662
663 proc_array = hypre_TAlloc(int, nprocs , HYPRE_MEMORY_HOST);
664 proc_array2 = hypre_TAlloc(int, nprocs , HYPRE_MEMORY_HOST);
665 for ( i = 0; i < nprocs; i++ ) proc_array2[i] = 0;
666 proc_array2[mypid] = Nrows;
667 MPI_Allreduce(proc_array2,proc_array,nprocs,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
668 NrowsOffset = 0;
669 for (i = 0; i < mypid; i++) NrowsOffset += proc_array[i];
670 for (i = 1; i < nprocs; i++) proc_array[i] += proc_array[i-1];
671 hypre_TFree(proc_array2, HYPRE_MEMORY_HOST);
672
673 /* ---------------------------------------------------------------- */
674 /* compose the column index map (index_array,index_array2) */
675 /* ---------------------------------------------------------------- */
676
677 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
678 context->comm = MPI_COMM_WORLD;
679 context->Amat = mh_mat;
680 dble_array = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
681 for (i = Nrows; i < extNrows; i++) dble_array[i] = 0.0;
682 for (i = 0; i < Nrows; i++) dble_array[i] = 1.0 * ( i + NrowsOffset );
683 MH_ExchBdry(dble_array, context);
684 if ( extNrows-Nrows > 0 )
685 index_array = hypre_TAlloc(int, (extNrows-Nrows) , HYPRE_MEMORY_HOST);
686 else
687 index_array = NULL;
688 for (i = Nrows; i < extNrows; i++) index_array[i-Nrows] = dble_array[i];
689 if ( extNrows-Nrows > 0 )
690 index_array2 = hypre_TAlloc(int, (extNrows-Nrows) , HYPRE_MEMORY_HOST);
691 else
692 index_array2 = NULL;
693 for (i = 0; i < extNrows-Nrows; i++) index_array2[i] = i;
694 hypre_TFree(dble_array, HYPRE_MEMORY_HOST);
695 hypre_TFree(context, HYPRE_MEMORY_HOST);
696
697 /* ---------------------------------------------------------------- */
698 /* send the lengths of each row to remote processor */
699 /* at the end, additional row information should be given */
700 /* in total_recv_leng, recv_lengths, int_buf, dble_buf */
701 /* ---------------------------------------------------------------- */
702
703 HYPRE_LSI_DDICTGetRowLengths(mh_mat, total_recv_leng, recv_lengths);
704 HYPRE_LSI_DDICTGetOffProcRows(mh_mat, *total_recv_leng, *recv_lengths,
705 NrowsOffset,index_array,index_array2,int_buf, dble_buf);
706
707 hypre_TFree(proc_array, HYPRE_MEMORY_HOST);
708 HYPRE_LSI_qsort1a(index_array, index_array2, 0, extNrows-Nrows-1);
709 (*sindex_array) = index_array;
710 (*sindex_array2) = index_array2;
711 (*offset) = NrowsOffset;
712 return 0;
713 }
714
715 /*****************************************************************************/
716 /* function for doing ICT decomposition */
717 /*---------------------------------------------------------------------------*/
718
HYPRE_LSI_DDICTDecompose(HYPRE_LSI_DDICT * ict_ptr,MH_Matrix * Amat,int total_recv_leng,int * recv_lengths,int * ext_ja,double * ext_aa,int * map,int * map2,int Noffset)719 int HYPRE_LSI_DDICTDecompose(HYPRE_LSI_DDICT *ict_ptr,MH_Matrix *Amat,
720 int total_recv_leng, int *recv_lengths, int *ext_ja, double *ext_aa,
721 int *map, int *map2, int Noffset)
722 {
723 int i, j, row_leng, *mat_ia, *mat_ja, allocated_space, *cols, mypid;
724 int index, ind2, total_nnz, offset, Nrows, extNrows;
725 double *vals, *mat_aa, *rowNorms, tau, rel_tau;
726 MH_Context *context;
727
728 /* ---------------------------------------------------------------- */
729 /* fetch ICT parameters */
730 /* ---------------------------------------------------------------- */
731
732 MPI_Comm_rank(ict_ptr->comm, &mypid);
733 tau = ict_ptr->thresh;
734 Nrows = Amat->Nrows;
735 extNrows = Nrows + total_recv_leng;
736 ict_ptr->Nrows = Nrows;
737 ict_ptr->extNrows = extNrows;
738
739 /* ---------------------------------------------------------------- */
740 /* allocate temporary storage space */
741 /* ---------------------------------------------------------------- */
742
743 allocated_space = extNrows;
744 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
745 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
746 rowNorms = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
747
748 /* ---------------------------------------------------------------- */
749 /* compute the storage requirement for the ILU matrix */
750 /* ---------------------------------------------------------------- */
751
752 context = hypre_TAlloc(MH_Context, 1, HYPRE_MEMORY_HOST);
753 context->Amat = Amat;
754 total_nnz = 0;
755 for ( i = 0; i < Nrows; i++ )
756 {
757 rowNorms[i] = 0.0;
758 while (MH_GetRow(context,1,&i,allocated_space,cols,vals,&row_leng)==0)
759 {
760 hypre_TFree(vals, HYPRE_MEMORY_HOST);
761 hypre_TFree(cols, HYPRE_MEMORY_HOST);
762 allocated_space += 200 + 1;
763 cols = hypre_TAlloc(int, allocated_space , HYPRE_MEMORY_HOST);
764 vals = hypre_TAlloc(double, allocated_space , HYPRE_MEMORY_HOST);
765 }
766 total_nnz += row_leng;
767 for ( j = 0; j < row_leng; j++ ) rowNorms[i] += habs(vals[j]);
768 rowNorms[i] /= extNrows;
769 rowNorms[i] = 1.0;
770 }
771 for ( i = 0; i < total_recv_leng; i++ ) total_nnz += recv_lengths[i];
772 mat_ia = hypre_TAlloc(int, (extNrows + 1 ) , HYPRE_MEMORY_HOST);
773 mat_ja = hypre_TAlloc(int, total_nnz , HYPRE_MEMORY_HOST);
774 mat_aa = hypre_TAlloc(double, total_nnz , HYPRE_MEMORY_HOST);
775
776 /* ---------------------------------------------------------------- */
777 /* construct the orginal matrix in CSR format */
778 /* ---------------------------------------------------------------- */
779
780 total_nnz = 0;
781 mat_ia[0] = 0;
782 for ( i = 0; i < Nrows; i++ )
783 {
784 rel_tau = tau * rowNorms[i];
785 MH_GetRow(context,1,&i,allocated_space,cols,vals,&row_leng);
786 for ( j = 0; j < row_leng; j++ )
787 {
788 if ( cols[j] <= i && habs(vals[j]) > rel_tau )
789 {
790 mat_aa[total_nnz] = vals[j];
791 mat_ja[total_nnz++] = cols[j];
792 }
793 }
794 mat_ia[i+1] = total_nnz;
795 }
796 offset = 0;
797 for ( i = 0; i < total_recv_leng; i++ )
798 {
799 rowNorms[i+Nrows] = 0.0;
800 for ( j = offset; j < offset+recv_lengths[i]; j++ )
801 {
802 index = ext_ja[j];
803 if ( index >= Noffset && index < Noffset+Nrows )
804 ext_ja[j] = index - Noffset;
805 else
806 {
807 ind2 = HYPRE_LSI_Search(map, index, extNrows-Nrows);
808 if ( ind2 >= 0 ) ext_ja[j] = map2[ind2] + Nrows;
809 else ext_ja[j] = -1;
810 }
811 if ( ext_ja[j] != -1 ) rowNorms[i+Nrows] += habs(ext_aa[j]);
812 }
813 rowNorms[i+Nrows] /= extNrows;
814 rowNorms[i+Nrows] = 1.0;
815 rel_tau = tau * rowNorms[i+Nrows];
816 for ( j = offset; j < offset+recv_lengths[i]; j++ )
817 {
818 if (ext_ja[j] != -1 && ext_ja[j] <= Nrows+i && habs(ext_aa[j]) > rel_tau)
819 {
820 mat_aa[total_nnz] = ext_aa[j];
821 mat_ja[total_nnz++] = ext_ja[j];
822 }
823 }
824 offset += recv_lengths[i];
825 mat_ia[Nrows+i+1] = total_nnz;
826 }
827
828 /* ---------------------------------------------------------------- */
829 /* clean up a little */
830 /* ---------------------------------------------------------------- */
831
832 hypre_TFree(Amat->rowptr, HYPRE_MEMORY_HOST);
833 hypre_TFree(Amat->colnum, HYPRE_MEMORY_HOST);
834 hypre_TFree(Amat->values, HYPRE_MEMORY_HOST);
835 hypre_TFree(context, HYPRE_MEMORY_HOST);
836 hypre_TFree(cols, HYPRE_MEMORY_HOST);
837 hypre_TFree(vals, HYPRE_MEMORY_HOST);
838
839 /* ---------------------------------------------------------------- */
840 /* call ICT factorization */
841 /* ---------------------------------------------------------------- */
842
843 HYPRE_LSI_DDICTFactorize(ict_ptr, mat_aa, mat_ja, mat_ia, rowNorms);
844
845 hypre_TFree(mat_aa , HYPRE_MEMORY_HOST);
846 hypre_TFree(mat_ia , HYPRE_MEMORY_HOST);
847 hypre_TFree(mat_ja , HYPRE_MEMORY_HOST);
848 hypre_TFree(rowNorms, HYPRE_MEMORY_HOST);
849
850 if ( ict_ptr->outputLevel > 0 )
851 {
852 total_nnz = ict_ptr->mat_ja[extNrows];
853 printf("%d : DDICT number of nonzeros = %d\n",mypid,total_nnz);
854 }
855
856 return 0;
857 }
858
859 /*****************************************************************************/
860 /* function for doing ICT factorization */
861 /*---------------------------------------------------------------------------*/
862
HYPRE_LSI_DDICTFactorize(HYPRE_LSI_DDICT * ict_ptr,double * mat_aa,int * mat_ja,int * mat_ia,double * rowNorms)863 int HYPRE_LSI_DDICTFactorize(HYPRE_LSI_DDICT *ict_ptr, double *mat_aa,
864 int *mat_ja, int *mat_ia, double *rowNorms)
865 {
866 int i, j, row_leng, first, row_beg, row_endp1, track_leng, *track_array;
867 int k, mypid, nnz_count, num_small_pivot, printstep, extNrows;
868 int *msr_iptr, *msc_jptr, *msc_jend, rowMax, Lcount, sortcnt, *sortcols;
869 int totalFill, colIndex, index;
870 double fillin, tau, rel_tau, *dble_buf, *msr_aptr, *msc_aptr, absval;
871 double *sortvals, ddata;
872
873 /* ---------------------------------------------------------------- */
874 /* fetch ICT parameters */
875 /* ---------------------------------------------------------------- */
876
877 MPI_Comm_rank(ict_ptr->comm, &mypid);
878 tau = ict_ptr->thresh;
879 fillin = ict_ptr->fillin;
880 extNrows = ict_ptr->extNrows;
881 rowMax = 0;
882 for ( i = 0; i < extNrows; i++ )
883 {
884 row_leng = mat_ia[i+1] - mat_ia[i];
885 if ( row_leng > rowMax ) rowMax = row_leng;
886 }
887 totalFill = rowMax * (fillin + 1) * extNrows;
888
889 /* ---------------------------------------------------------------- */
890 /* allocate permanent and temporary storage */
891 /* ---------------------------------------------------------------- */
892
893 track_array = hypre_TAlloc(int, extNrows , HYPRE_MEMORY_HOST);
894 sortcols = hypre_TAlloc(int, extNrows , HYPRE_MEMORY_HOST);
895 sortvals = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
896 dble_buf = hypre_TAlloc(double, extNrows , HYPRE_MEMORY_HOST);
897 msr_iptr = hypre_TAlloc(int, (totalFill+extNrows+1) , HYPRE_MEMORY_HOST);
898 msc_jptr = hypre_TAlloc(int, (totalFill+extNrows+1) , HYPRE_MEMORY_HOST);
899 msc_jend = hypre_TAlloc(int, (extNrows + 1 ) , HYPRE_MEMORY_HOST);
900 msr_aptr = hypre_TAlloc(double, (totalFill+extNrows) , HYPRE_MEMORY_HOST);
901 msc_aptr = hypre_TAlloc(double, (totalFill+extNrows) , HYPRE_MEMORY_HOST);
902 msc_jptr[0] = msc_jend[0] = extNrows + 1;
903 for ( i = 1; i <= extNrows; i++ )
904 {
905 msc_jptr[i] = msc_jptr[i-1] + rowMax * (fillin + 1);
906 msc_jend[i] = msc_jptr[i];
907 }
908 for ( i = 0; i < extNrows; i++ ) dble_buf[i] = 0.0;
909 printstep = extNrows / 10;
910
911 /* ---------------------------------------------------------------- */
912 /* process the rows */
913 /* ---------------------------------------------------------------- */
914
915 num_small_pivot = 0;
916 nnz_count = extNrows + 1;
917 msr_iptr[0] = nnz_count;
918
919 for ( i = 0; i < extNrows; i++ )
920 {
921 if ( i % printstep == 0 && ict_ptr->outputLevel > 0 )
922 printf("%4d : DDICT Processing row %6d (%6d)\n",mypid,i,extNrows);
923
924 /* ------------------------------------------------------------- */
925 /* get the row information */
926 /* ------------------------------------------------------------- */
927
928 track_leng = 0;
929 row_beg = mat_ia[i];
930 row_endp1 = mat_ia[i+1];
931 row_leng = row_endp1 - row_beg;
932 first = i;
933 rel_tau = tau * rowNorms[i];
934
935 /* ------------------------------------------------------------- */
936 /* load the row into dble_buf */
937 /* ------------------------------------------------------------- */
938
939 for ( j = row_beg; j < row_endp1; j++ )
940 {
941 colIndex = mat_ja[j];
942 if ( colIndex > i ) printf("WARNING (A)\n");
943 dble_buf[colIndex] = mat_aa[j];
944 track_array[track_leng++] = colIndex;
945 if ( colIndex < first ) first = colIndex;
946 }
947 Lcount = row_leng * fillin;
948
949 /* ------------------------------------------------------------- */
950 /* reduce the row */
951 /* ------------------------------------------------------------- */
952
953 for ( j = first; j < i; j++ )
954 {
955 if ( habs(dble_buf[j]) > rel_tau )
956 {
957 ddata = dble_buf[j] * msr_aptr[j];
958
959 for ( k = msc_jptr[j]; k < msc_jend[j]; k++ )
960 {
961 colIndex = msc_jptr[k];
962 if ( colIndex > j && colIndex < i )
963 {
964 if ( dble_buf[colIndex] != 0.0 )
965 dble_buf[colIndex] -= (ddata * msc_aptr[k]);
966 else
967 {
968 dble_buf[colIndex] = - (ddata * msc_aptr[k]);
969 track_array[track_leng++] = colIndex;
970 }
971 }
972 }
973 dble_buf[j] = ddata;
974 }
975 else dble_buf[j] = 0.0;
976 }
977
978 /* ------------------------------------------------------------- */
979 /* sort the new nonzeros */
980 /* ------------------------------------------------------------- */
981
982 sortcnt = 0;
983 if ( track_leng > extNrows ) printf("WARNING (B)\n");
984 for ( j = row_leng; j < track_leng; j++ )
985 {
986 index = track_array[j];
987 absval = habs(dble_buf[index]);
988 if ( absval > rel_tau )
989 {
990 sortcols[sortcnt] = index;
991 sortvals[sortcnt++] = absval * rowNorms[index];
992 }
993 else dble_buf[index] = 0.0;
994 }
995 if ( sortcnt > Lcount )
996 {
997 HYPRE_LSI_SplitDSort(sortvals,sortcnt,sortcols,Lcount);
998 for ( j = Lcount; j < sortcnt; j++ ) dble_buf[sortcols[j]] = 0.0;
999 for ( j = 0; j < row_leng; j++ )
1000 {
1001 index = track_array[j];
1002 if ( index != i )
1003 {
1004 ddata = dble_buf[i] - (dble_buf[index] * dble_buf[index]);
1005 if ( ddata > 1.0E-10 ) dble_buf[i] = ddata;
1006 else
1007 {
1008 printf("%d : DDICT negative pivot (%d,%d,%d)\n", mypid,
1009 i, j, extNrows);
1010 num_small_pivot++;
1011 for ( k = j; k < row_leng; k++ )
1012 {
1013 index = track_array[k];
1014 dble_buf[index] = 0.0;
1015 }
1016 Lcount = 0;
1017 break;
1018 }
1019 }
1020 }
1021 for ( j = 0; j < Lcount; j++ )
1022 {
1023 index = sortcols[j];
1024 ddata = dble_buf[i] - (dble_buf[index] * dble_buf[index]);
1025 if ( ddata > 1.0E-10 ) dble_buf[i] = ddata;
1026 else
1027 {
1028 printf("%d : (2) DDICT negative pivot (%d,%d,%d)\n", mypid,
1029 i, j, extNrows);
1030 num_small_pivot++;
1031 for ( k = j; k < Lcount; k++ )
1032 {
1033 index = sortcols[k];
1034 dble_buf[index] = 0.0;
1035 }
1036 Lcount = j;
1037 break;
1038 }
1039 }
1040 }
1041 else
1042 {
1043 for ( j = 0; j < row_leng; j++ )
1044 {
1045 index = track_array[j];
1046 if ( index != i )
1047 {
1048 ddata = dble_buf[i] - (dble_buf[index] * dble_buf[index]);
1049 if ( ddata > 1.0E-10 ) dble_buf[i] = ddata;
1050 else
1051 {
1052 printf("%d : DDICT negative pivot (%d,%d,%d)\n", mypid,
1053 i, j, extNrows);
1054 num_small_pivot++;
1055 for ( k = j; k < row_leng; k++ )
1056 {
1057 index = track_array[k];
1058 dble_buf[index] = 0.0;
1059 }
1060 sortcnt = 0;
1061 break;
1062 }
1063 }
1064 }
1065 for ( j = 0; j < sortcnt; j++ )
1066 {
1067 index = sortcols[j];
1068 ddata = dble_buf[i] - (dble_buf[index] * dble_buf[index]);
1069 if ( ddata > 1.0E-10 ) dble_buf[i] = ddata;
1070 else
1071 {
1072 printf("%d : (2) DDICT negative pivot (%d,%d,%d)\n", mypid,
1073 i, j, extNrows);
1074 num_small_pivot++;
1075 for ( k = j; k < sortcnt; k++ )
1076 {
1077 index = sortcols[k];
1078 dble_buf[index] = 0.0;
1079 }
1080 sortcnt = j;
1081 break;
1082 }
1083 }
1084 }
1085 if ( dble_buf[i] > 0 )
1086 {
1087 if ( dble_buf[i] < 1.0E-10 )
1088 {
1089 num_small_pivot++;
1090 msc_aptr[i] = msr_aptr[i] = 1.0E5;
1091 }
1092 else msc_aptr[i] = msr_aptr[i] = 1.0 / sqrt( dble_buf[i] );
1093 dble_buf[i] = 0.0;
1094 }
1095 else
1096 {
1097 printf("%4d : ERROR in DDICT - negative or zero pivot.\n", mypid);
1098 printf(" L(%4d,%4d) = %e\n", i, i, dble_buf[i]);
1099 msc_aptr[i] = msr_aptr[i] = 1.0 / sqrt( - dble_buf[i] );
1100 dble_buf[i] = 0.0;
1101 }
1102 for ( j = 0; j < track_leng; j++ )
1103 {
1104 index = track_array[j];
1105 if ( index < i && dble_buf[index] != 0.0 )
1106 {
1107 msr_aptr[nnz_count] = dble_buf[index];
1108 msr_iptr[nnz_count++] = index;
1109 colIndex = msc_jend[index]++;
1110 msc_aptr[colIndex] = dble_buf[index];
1111 msc_jptr[colIndex] = i;
1112 dble_buf[index] = 0.0;
1113 }
1114 }
1115 msr_iptr[i+1] = nnz_count;
1116 }
1117
1118 if ( nnz_count > totalFill+extNrows )
1119 printf("%4d : DDICT WARNING : buffer overflow (%d,%d)\n",mypid,nnz_count,
1120 totalFill+extNrows);
1121 if ( ict_ptr->outputLevel > 0 )
1122 {
1123 printf("%4d : DDICT number of nonzeros = %d\n",mypid,nnz_count);
1124 printf("%4d : DDICT number of small pivots = %d\n",mypid,num_small_pivot);
1125 }
1126
1127 /* ---------------------------------------------------------- */
1128 /* deallocate temporary storage space */
1129 /* ---------------------------------------------------------- */
1130
1131 hypre_TFree(track_array, HYPRE_MEMORY_HOST);
1132 hypre_TFree(sortcols, HYPRE_MEMORY_HOST);
1133 hypre_TFree(sortvals, HYPRE_MEMORY_HOST);
1134 hypre_TFree(dble_buf, HYPRE_MEMORY_HOST);
1135 hypre_TFree(msc_jptr, HYPRE_MEMORY_HOST);
1136 hypre_TFree(msc_jend, HYPRE_MEMORY_HOST);
1137 hypre_TFree(msc_aptr, HYPRE_MEMORY_HOST);
1138
1139 ict_ptr->mat_ja = msr_iptr;
1140 ict_ptr->mat_aa = msr_aptr;
1141 return 0;
1142 }
1143
1144