43 #include "Factor_dh.h"
46 #include "SubdomainGraph_dh.h"
47 #include "TimeLog_dh.h"
49 #include "Numbering_dh.h"
50 #include "Hash_i_dh.h"
51 #include "Parser_dh.h"
52 #include "mat_dh_private.h"
53 #include "getRow_dh.h"
54 #include "Euclid_dh.h"
63 static void adjust_bj_private (
Factor_dh mat);
64 static void unadjust_bj_private (
Factor_dh mat);
68 #define __FUNC__ "Factor_dhCreate"
74 if (np_dh > MAX_MPI_TASKS)
76 SET_V_ERROR (
"you must change MAX_MPI_TASKS and recompile!");
89 tmp->blockJacobi =
false;
98 tmp->work_y_lo = tmp->work_x_hi = NULL;
99 tmp->sendbufLo = tmp->sendbufHi = NULL;
100 tmp->sendindLo = tmp->sendindHi = NULL;
101 tmp->num_recvLo = tmp->num_recvHi = 0;
102 tmp->num_sendLo = tmp->num_sendHi = 0;
103 tmp->sendlenLo = tmp->sendlenHi = 0;
105 tmp->solveIsSetup =
false;
106 tmp->numbSolve = NULL;
108 tmp->debug = Parser_dhHasSwitch (parser_dh,
"-debug_Factor");
114 #define __FUNC__ "Factor_dhDestroy"
118 START_FUNC_DH
if (mat->rp != NULL)
123 if (mat->cval != NULL)
128 if (mat->aval != NULL)
133 if (mat->diag != NULL)
138 if (mat->fill != NULL)
144 if (mat->work_y_lo != NULL)
146 FREE_DH (mat->work_y_lo);
149 if (mat->work_x_hi != NULL)
151 FREE_DH (mat->work_x_hi);
154 if (mat->sendbufLo != NULL)
156 FREE_DH (mat->sendbufLo);
159 if (mat->sendbufHi != NULL)
161 FREE_DH (mat->sendbufHi);
164 if (mat->sendindLo != NULL)
166 FREE_DH (mat->sendindLo);
169 if (mat->sendindHi != NULL)
171 FREE_DH (mat->sendindHi);
175 if (mat->numbSolve != NULL)
177 Numbering_dhDestroy (mat->numbSolve);
186 #define __FUNC__ "create_fake_mat_private"
190 START_FUNC_DH
Mat_dh matFake;
191 Mat_dhCreate (matFakeIN);
193 matFake = *matFakeIN;
196 matFake->rp = mat->rp;
197 matFake->cval = mat->cval;
198 matFake->aval = mat->aval;
199 matFake->beg_row = mat->beg_row;
203 #define __FUNC__ "destroy_fake_mat_private"
205 destroy_fake_mat_private (
Mat_dh matFake)
207 START_FUNC_DH matFake->rp = NULL;
208 matFake->cval = NULL;
209 matFake->aval = NULL;
210 Mat_dhDestroy (matFake);
217 #define __FUNC__ "Factor_dhReadNz"
221 START_FUNC_DH
int ierr, retval = mat->rp[mat->m];
223 ierr = MPI_Allreduce (&nz, &retval, 1, MPI_INT, MPI_SUM, comm_dh);
224 CHECK_MPI_ERROR (ierr);
225 END_FUNC_VAL (retval)}
230 #define __FUNC__ "Factor_dhPrintRows"
232 Factor_dhPrintRows (
Factor_dh mat, FILE * fp)
234 START_FUNC_DH
int beg_row = mat->beg_row;
235 int m = mat->m, i, j;
238 noValues = (Parser_dhHasSwitch (parser_dh,
"-noValues"));
239 if (mat->aval == NULL)
242 if (mat->blockJacobi)
244 adjust_bj_private (mat);
249 "\n----------------------- Factor_dhPrintRows ------------------\n");
250 if (mat->blockJacobi)
253 "@@@ Block Jacobi ILU; adjusted values from zero-based @@@\n");
256 for (i = 0; i < m; ++i)
258 fprintf (fp,
"%i :: ", 1 + i + beg_row);
259 for (j = mat->rp[i]; j < mat->rp[i + 1]; ++j)
263 fprintf (fp,
"%i ", 1 + mat->cval[j]);
267 fprintf (fp,
"%i,%g ; ", 1 + mat->cval[j], mat->aval[j]);
273 if (mat->blockJacobi)
275 unadjust_bj_private (mat);
282 #define __FUNC__ "Factor_dhPrintDiags"
284 Factor_dhPrintDiags (
Factor_dh mat, FILE * fp)
286 START_FUNC_DH
int beg_row = mat->beg_row;
287 int m = mat->m, i, pe, *diag = mat->diag;
288 REAL_DH *aval = mat->aval;
292 "\n----------------------- Factor_dhPrintDiags ------------------\n");
293 fprintf_dh (fp,
"(grep for 'ZERO')\n");
295 for (pe = 0; pe < np_dh; ++pe)
297 MPI_Barrier (comm_dh);
300 fprintf (fp,
"----- subdomain: %i processor: %i\n", pe, myid_dh);
301 for (i = 0; i < m; ++i)
303 REAL_DH val = aval[diag[i]];
306 fprintf (fp,
"%i %g\n", i + 1 + beg_row, aval[diag[i]]);
310 fprintf (fp,
"%i %g ZERO\n", i + 1 + beg_row,
320 #define __FUNC__ "Factor_dhPrintGraph"
322 Factor_dhPrintGraph (
Factor_dh mat,
char *filename)
324 START_FUNC_DH FILE * fp;
325 int i, j, m = mat->m, *work, *rp = mat->rp, *cval = mat->cval;
328 SET_V_ERROR (
"only implemented for single mpi task");
330 work = (
int *) MALLOC_DH (m *
sizeof (
int));
333 fp = openFile_dh (filename,
"w");
336 for (i = 0; i < m; ++i)
338 for (j = 0; j < m; ++j)
340 for (j = rp[i]; j < rp[i]; ++j)
343 for (j = 0; j < m; ++j)
365 #define __FUNC__ "Factor_dhPrintTriples"
367 Factor_dhPrintTriples (
Factor_dh mat,
char *filename)
369 START_FUNC_DH
int pe, i, j;
370 int m = mat->m, *rp = mat->rp;
371 int beg_row = mat->beg_row;
372 REAL_DH *aval = mat->aval;
376 if (mat->blockJacobi)
378 adjust_bj_private (mat);
382 noValues = (Parser_dhHasSwitch (parser_dh,
"-noValues"));
386 for (pe = 0; pe < np_dh; ++pe)
388 MPI_Barrier (comm_dh);
393 fp = openFile_dh (filename,
"w");
398 fp = openFile_dh (filename,
"a");
402 for (i = 0; i < m; ++i)
404 for (j = rp[i]; j < rp[i + 1]; ++j)
408 fprintf (fp,
"%i %i\n", 1 + i + beg_row,
413 fprintf (fp, TRIPLES_FORMAT,
414 1 + i + beg_row, 1 + mat->cval[j], aval[j]);
423 if (mat->blockJacobi)
425 unadjust_bj_private (mat);
446 #define __FUNC__ "setup_receives_private"
448 setup_receives_private (
Factor_dh mat,
int *beg_rows,
int *end_rows,
449 double *recvBuf, MPI_Request * req,
450 int *reqind,
int reqlen,
int *outlist,
bool debug)
452 START_FUNC_DH
int i, j, this_pe, num_recv = 0;
458 "\nFACT ========================================================\n");
459 fprintf (logFile,
"FACT STARTING: setup_receives_private\n");
462 for (i = 0; i < reqlen; i = j)
465 this_pe = mat_find_owner (beg_rows, end_rows, reqind[i]);
469 for (j = i + 1; j < reqlen; j++)
472 if (idx < beg_rows[this_pe] || idx >= end_rows[this_pe])
481 fprintf (logFile,
"FACT need nodes from P_%i: ", this_pe);
482 for (k = i; k < j; ++k)
483 fprintf (logFile,
"%i ", 1 + reqind[k]);
484 fprintf (logFile,
"\n");
488 outlist[this_pe] = j - i;
496 MPI_Isend (reqind + i, j - i, MPI_INT, this_pe, 444, comm_dh, &request);
497 MPI_Request_free (&request);
500 MPI_Recv_init (recvBuf + i, j - i, MPI_DOUBLE, this_pe, 555,
501 comm_dh, req + num_recv);
505 END_FUNC_VAL (num_recv);
514 #define __FUNC__ "setup_sends_private"
516 setup_sends_private (
Factor_dh mat,
int *inlist,
517 int *o2n_subdomain,
bool debug)
519 START_FUNC_DH
int i, jLo, jHi, sendlenLo, sendlenHi, first = mat->beg_row;
520 MPI_Request *requests = mat->requests, *sendReq;
521 MPI_Status *statuses = mat->status;
525 int myidNEW = o2n_subdomain[myid_dh];
530 fprintf (logFile,
"FACT \nSTARTING: setup_sends_private\n");
534 sendlenLo = sendlenHi = 0;
535 for (i = 0; i < np_dh; i++)
539 if (o2n_subdomain[i] < myidNEW)
541 sendlenLo += inlist[i];
545 sendlenHi += inlist[i];
550 mat->sendlenLo = sendlenLo;
551 mat->sendlenHi = sendlenHi;
552 mat->sendbufLo = (
double *) MALLOC_DH (sendlenLo *
sizeof (
double));
554 mat->sendbufHi = (
double *) MALLOC_DH (sendlenHi *
sizeof (
double));
556 mat->sendindLo = (
int *) MALLOC_DH (sendlenLo *
sizeof (
int));
558 mat->sendindHi = (
int *) MALLOC_DH (sendlenHi *
sizeof (
int));
565 for (i = 0; i < np_dh; i++)
569 isHigher = (o2n_subdomain[i] < myidNEW) ?
false :
true;
574 rcvBuf = &mat->sendindHi[jHi];
575 sendBuf = &mat->sendbufHi[jHi];
576 sendReq = &mat->send_reqHi[mat->num_sendHi];
582 rcvBuf = &mat->sendindLo[jLo];
583 sendBuf = &mat->sendbufLo[jLo];
584 sendReq = &mat->send_reqLo[mat->num_sendLo];
592 MPI_Irecv (rcvBuf, inlist[i], MPI_INT, i, 444, comm_dh,
597 MPI_Send_init (sendBuf, inlist[i], MPI_DOUBLE, i, 555, comm_dh,
603 MPI_Waitall (count, requests, statuses);
611 "\nFACT columns that I must send to other subdomains:\n");
612 for (i = 0; i < np_dh; i++)
616 isHigher = (o2n_subdomain[i] < myidNEW) ?
false :
true;
619 rcvBuf = &mat->sendindHi[jHi];
624 rcvBuf = &mat->sendindLo[jLo];
628 fprintf (logFile,
"FACT send to P_%i: ", i);
629 for (j = 0; j < inlist[i]; ++j)
630 fprintf (logFile,
"%i ", rcvBuf[j] + 1);
631 fprintf (logFile,
"\n");
638 for (i = 0; i < mat->sendlenLo; i++)
639 mat->sendindLo[i] -= first;
640 for (i = 0; i < mat->sendlenHi; i++)
641 mat->sendindHi[i] -= first;
647 #define __FUNC__ "Factor_dhSolveSetup"
651 START_FUNC_DH
int *outlist, *inlist;
652 int i, row, *rp = mat->rp, *cval = mat->cval;
656 int *beg_rows = sg->beg_rowP, *row_count = sg->row_count, *end_rows;
661 if (mat->debug && logFile != NULL)
664 end_rows = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
666 outlist = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
668 inlist = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
670 for (i = 0; i < np_dh; ++i)
674 end_rows[i] = beg_rows[i] + row_count[i];
678 create_fake_mat_private (mat, &matFake);
680 Numbering_dhCreate (&(mat->numbSolve));
682 numb = mat->numbSolve;
683 Numbering_dhSetup (numb, matFake);
685 destroy_fake_mat_private (matFake);
690 fprintf (stderr,
"Numbering_dhSetup completed\n");
694 i = m + numb->num_ext;
695 mat->work_y_lo = (
double *) MALLOC_DH (i *
sizeof (
double));
697 mat->work_x_hi = (
double *) MALLOC_DH (i *
sizeof (
double));
701 fprintf (logFile,
"FACT num_extLo= %i num_extHi= %i\n",
702 numb->num_extLo, numb->num_extHi);
709 recvBuf = mat->work_y_lo + m;
710 mat->num_recvLo = setup_receives_private (mat, beg_rows, end_rows,
711 recvBuf, mat->recv_reqLo,
713 numb->num_extLo, outlist,
721 recvBuf = mat->work_x_hi + m + numb->num_extLo;
722 mat->num_recvHi = setup_receives_private (mat, beg_rows, end_rows,
723 recvBuf, mat->recv_reqHi,
725 numb->num_extHi, outlist,
730 MPI_Alltoall (outlist, 1, MPI_INT, inlist, 1, MPI_INT, comm_dh);
737 setup_sends_private (mat, inlist, sg->o2n_sub, debug);
741 for (row = 0; row < m; row++)
743 int len = rp[row + 1] - rp[row];
744 int *ind = cval + rp[row];
745 Numbering_dhGlobalToLocal (numb, len, ind, ind);
761 "\n--------- row/col structure, after global to local renumbering\n");
762 for (ii = 0; ii < mat->m; ++ii)
764 fprintf (logFile,
"local row %i :: ", ii + 1);
765 for (jj = mat->rp[ii]; jj < mat->rp[ii + 1]; ++jj)
767 fprintf (logFile,
"%i ", 1 + mat->cval[jj]);
769 fprintf (logFile,
"\n");
771 fprintf (logFile,
"\n");
780 static void forward_solve_private (
int m,
int from,
int to,
781 int *rp,
int *cval,
int *diag,
782 double *aval,
double *rhs,
double *work_y,
785 static void backward_solve_private (
int m,
int from,
int to,
786 int *rp,
int *cval,
int *diag,
787 double *aval,
double *work_y,
788 double *work_x,
bool debug);
794 #define __FUNC__ "Factor_dhSolve"
796 Factor_dhSolve (
double *rhs,
double *lhs,
Euclid_dh ctx)
800 int ierr, i, m = mat->m, first_bdry = mat->first_bdry;
801 int offsetLo = mat->numbSolve->num_extLo;
802 int offsetHi = mat->numbSolve->num_extHi;
803 int *rp = mat->rp, *cval = mat->cval, *diag = mat->diag;
804 double *aval = mat->aval;
805 int *sendindLo = mat->sendindLo, *sendindHi = mat->sendindHi;
806 int sendlenLo = mat->sendlenLo, sendlenHi = mat->sendlenHi;
807 double *sendbufLo = mat->sendbufLo, *sendbufHi = mat->sendbufHi;
808 double *work_y = mat->work_y_lo;
809 double *work_x = mat->work_x_hi;
812 if (mat->debug && logFile != NULL)
815 beg_rowG = ctx->F->beg_row;
827 "\n=====================================================\n");
829 "FACT Factor_dhSolve: num_recvLo= %i num_recvHi = %i\n",
830 mat->num_recvLo, mat->num_recvHi);
836 MPI_Startall (mat->num_recvLo, mat->recv_reqLo);
840 MPI_Startall (mat->num_recvHi, mat->recv_reqHi);
851 forward_solve_private (m, from, to, rp, cval, diag, aval,
861 MPI_Waitall (mat->num_recvLo, mat->recv_reqLo, mat->status);
867 "FACT got 'y' values from lower neighbors; work buffer:\n ");
868 for (i = 0; i < offsetLo; ++i)
870 fprintf (logFile,
"%g ", work_y[m + i]);
880 forward_solve_private (m, from, to, rp, cval, diag, aval,
890 for (i = 0; i < sendlenHi; i++)
892 sendbufHi[i] = work_y[sendindHi[i]];
896 MPI_Startall (mat->num_sendHi, mat->send_reqHi);
902 "\nFACT sending 'y' values to higher neighbor:\nFACT ");
903 for (i = 0; i < sendlenHi; i++)
905 fprintf (logFile,
"%g ", sendbufHi[i]);
907 fprintf (logFile,
"\n");
917 ierr = MPI_Waitall (mat->num_recvHi, mat->recv_reqHi, mat->status);
918 CHECK_MPI_V_ERROR (ierr);
923 fprintf (logFile,
"FACT got 'x' values from higher neighbors:\n ");
924 for (i = m + offsetLo; i < m + offsetLo + offsetHi; ++i)
926 fprintf (logFile,
"%g ", work_x[i]);
928 fprintf (logFile,
"\n");
937 backward_solve_private (m, from, to, rp, cval, diag, aval,
938 work_y, work_x, debug);
947 for (i = 0; i < sendlenLo; i++)
949 sendbufLo[i] = work_x[sendindLo[i]];
953 ierr = MPI_Startall (mat->num_sendLo, mat->send_reqLo);
954 CHECK_MPI_V_ERROR (ierr);
960 "\nFACT sending 'x' values to lower neighbor:\nFACT ");
961 for (i = 0; i < sendlenLo; i++)
963 fprintf (logFile,
"%g ", sendbufLo[i]);
965 fprintf (logFile,
"\n");
974 backward_solve_private (m, from, to, rp, cval, diag, aval,
975 work_y, work_x, debug);
980 memcpy (lhs, work_x, m *
sizeof (
double));
984 fprintf (logFile,
"\nFACT solution: ");
985 for (i = 0; i < m; ++i)
987 fprintf (logFile,
"%g ", lhs[i]);
989 fprintf (logFile,
"\n");
995 ierr = MPI_Waitall (mat->num_sendLo, mat->send_reqLo, mat->status);
996 CHECK_MPI_V_ERROR (ierr);
1001 ierr = MPI_Waitall (mat->num_sendHi, mat->send_reqHi, mat->status);
1002 CHECK_MPI_V_ERROR (ierr);
1009 #define __FUNC__ "forward_solve_private"
1011 forward_solve_private (
int m,
int from,
int to,
int *rp,
1012 int *cval,
int *diag,
double *aval,
1013 double *rhs,
double *work_y,
bool debug)
1015 START_FUNC_DH
int i, j, idx;
1020 "\nFACT starting forward_solve_private; from= %i; to= %i, m= %i\n",
1021 1 + from, 1 + to, m);
1037 for (i = from; i < to; ++i)
1039 int len = diag[i] - rp[i];
1040 int *col = cval + rp[i];
1041 double *val = aval + rp[i];
1042 double sum = rhs[i];
1044 fprintf (logFile,
"FACT solving for work_y[%i] (global)\n",
1046 fprintf (logFile,
"FACT sum = %g\n", sum);
1047 for (j = 0; j < len; ++j)
1050 sum -= (val[j] * work_y[idx]);
1052 "FACT sum(%g) -= val[j] (%g) * work_y[%i] (%g)\n",
1053 sum, val[j], 1 + idx, work_y[idx]);
1056 fprintf (logFile,
"FACT work_y[%i] = %g\n", 1 + i + beg_rowG,
1058 fprintf (logFile,
"-----------\n");
1061 fprintf (logFile,
"\nFACT work vector at end of forward solve:\n");
1062 for (i = 0; i < to; i++)
1063 fprintf (logFile,
" %i %g\n", i + 1 + beg_rowG, work_y[i]);
1068 for (i = from; i < to; ++i)
1070 int len = diag[i] - rp[i];
1071 int *col = cval + rp[i];
1072 double *val = aval + rp[i];
1073 double sum = rhs[i];
1075 for (j = 0; j < len; ++j)
1078 sum -= (val[j] * work_y[idx]);
1086 #define __FUNC__ "backward_solve_private"
1088 backward_solve_private (
int m,
int from,
int to,
int *rp,
1089 int *cval,
int *diag,
double *aval,
1090 double *work_y,
double *work_x,
bool debug)
1092 START_FUNC_DH
int i, j, idx;
1097 "\nFACT starting backward_solve_private; from= %i; to= %i, m= %i\n",
1098 1 + from, 1 + to, m);
1099 for (i = from - 1; i >= to; --i)
1101 int len = rp[i + 1] - diag[i] - 1;
1102 int *col = cval + diag[i] + 1;
1103 double *val = aval + diag[i] + 1;
1104 double sum = work_y[i];
1105 fprintf (logFile,
"FACT solving for work_x[%i]\n",
1108 for (j = 0; j < len; ++j)
1111 sum -= (val[j] * work_x[idx]);
1113 "FACT sum(%g) -= val[j] (%g) * work_x[idx] (%g)\n",
1114 sum, val[j], work_x[idx]);
1116 work_x[i] = sum * aval[diag[i]];
1117 fprintf (logFile,
"FACT work_x[%i] = %g\n", 1 + i, work_x[i]);
1118 fprintf (logFile,
"----------\n");
1124 for (i = from - 1; i >= to; --i)
1126 int len = rp[i + 1] - diag[i] - 1;
1127 int *col = cval + diag[i] + 1;
1128 double *val = aval + diag[i] + 1;
1129 double sum = work_y[i];
1131 for (j = 0; j < len; ++j)
1134 sum -= (val[j] * work_x[idx]);
1136 work_x[i] = sum * aval[diag[i]];
1142 #define __FUNC__ "Factor_dhInit"
1144 Factor_dhInit (
void *A,
bool fillFlag,
bool avalFlag,
1145 double rho,
int id,
int beg_rowP,
Factor_dh * Fout)
1147 START_FUNC_DH
int m, n, beg_row, alloc;
1150 EuclidGetDimensions (A, &beg_row, &m, &n);
1153 Factor_dhCreate (&F);
1159 F->beg_row = beg_rowP;
1163 F->rp = (
int *) MALLOC_DH ((m + 1) *
sizeof (int));
1166 F->cval = (
int *) MALLOC_DH (alloc *
sizeof (
int));
1168 F->diag = (
int *) MALLOC_DH (m *
sizeof (
int));
1172 F->fill = (
int *) MALLOC_DH (alloc *
sizeof (
int));
1177 F->aval = (REAL_DH *) MALLOC_DH (alloc *
sizeof (REAL_DH));
1183 #define __FUNC__ "Factor_dhReallocate"
1185 Factor_dhReallocate (
Factor_dh F,
int used,
int additional)
1187 START_FUNC_DH
int alloc = F->alloc;
1189 if (used + additional > F->alloc)
1192 while (alloc < used + additional)
1196 F->cval = (
int *) MALLOC_DH (alloc *
sizeof (
int));
1198 memcpy (F->cval, tmpI, used * sizeof (
int));
1201 if (F->fill != NULL)
1204 F->fill = (
int *) MALLOC_DH (alloc *
sizeof (
int));
1206 memcpy (F->fill, tmpI, used * sizeof (
int));
1210 if (F->aval != NULL)
1212 REAL_DH *tmpF = F->aval;
1213 F->aval = (REAL_DH *) MALLOC_DH (alloc *
sizeof (REAL_DH));
1215 memcpy (F->aval, tmpF, used * sizeof (REAL_DH));
1223 #define __FUNC__ "Factor_dhTranspose"
1231 SET_V_ERROR (
"only for sequential");
1234 Factor_dhCreate (&B);
1238 if (B->aval == NULL)
1240 mat_dh_transpose_private (A->m, A->rp, &B->rp, A->cval, &B->cval,
1246 mat_dh_transpose_private (A->m, A->rp, &B->rp, A->cval, &B->cval,
1255 #define __FUNC__ "Factor_dhSolveSeq"
1257 Factor_dhSolveSeq (
double *rhs,
double *lhs,
Euclid_dh ctx)
1262 printf (
"F is null.\n");
1266 printf (
"F isn't null.\n");
1268 int *rp, *cval, *diag;
1269 int i, j, *vi, nz, m = F->m;
1270 REAL_DH *aval, *work;
1275 if (ctx->F->debug && logFile != NULL)
1288 "\nFACT ============================================================\n");
1289 fprintf (logFile,
"FACT starting Factor_dhSolveSeq\n");
1292 fprintf (logFile,
"\nFACT STARTING FORWARD SOLVE\n------------\n");
1294 fprintf (logFile,
"FACT work[0] = %g\n------------\n", work[0]);
1295 for (i = 1; i < m; i++)
1299 nz = diag[i] - rp[i];
1300 fprintf (logFile,
"FACT solving for work[%i]\n", i + 1);
1302 for (j = 0; j < nz; ++j)
1304 sum -= (v[j] * work[vi[j]]);
1306 "FACT sum (%g) -= v[j] (%g) * work[vi[j]] (%g)\n",
1307 sum, v[j], work[vi[j]]);
1310 fprintf (logFile,
"FACT work[%i] = %g\n------------\n", 1 + i,
1315 fprintf (logFile,
"\nFACT work vector at end of forward solve:\n");
1316 for (i = 0; i < m; i++)
1317 fprintf (logFile,
" %i %g\n", i + 1, work[i]);
1321 fprintf (logFile,
"\nFACT STARTING BACKWARD SOLVE\n--------------\n");
1322 for (i = m - 1; i >= 0; i--)
1324 v = aval + diag[i] + 1;
1325 vi = cval + diag[i] + 1;
1326 nz = rp[i + 1] - diag[i] - 1;
1327 fprintf (logFile,
"FACT solving for lhs[%i]\n", i + 1);
1329 for (j = 0; j < nz; ++j)
1331 sum -= (v[j] * work[vi[j]]);
1333 "FACT sum (%g) -= v[j] (%g) * work[vi[j]] (%g)\n",
1334 sum, v[j], work[vi[j]]);
1336 lhs[i] = work[i] = sum * aval[diag[i]];
1337 fprintf (logFile,
"FACT lhs[%i] = %g\n------------\n", 1 + i,
1339 fprintf (logFile,
"FACT solving for lhs[%i]\n", i + 1);
1342 fprintf (logFile,
"\nFACT solution: ");
1343 for (i = 0; i < m; ++i)
1344 fprintf (logFile,
"%g ", lhs[i]);
1345 fprintf (logFile,
"\n");
1353 for (i = 1; i < m; i++)
1357 nz = diag[i] - rp[i];
1360 sum -= (*v++ * work[*vi++]);
1365 for (i = m - 1; i >= 0; i--)
1367 v = aval + diag[i] + 1;
1368 vi = cval + diag[i] + 1;
1369 nz = rp[i + 1] - diag[i] - 1;
1372 sum -= (*v++ * work[*vi++]);
1373 lhs[i] = work[i] = sum * aval[diag[i]];
1383 #define __FUNC__ "adjust_bj_private"
1387 START_FUNC_DH
int i;
1388 int nz = mat->rp[mat->m];
1389 int beg_row = mat->beg_row;
1390 for (i = 0; i < nz; ++i)
1391 mat->cval[i] += beg_row;
1395 #define __FUNC__ "unadjust_bj_private"
1399 START_FUNC_DH
int i;
1400 int nz = mat->rp[mat->m];
1401 int beg_row = mat->beg_row;
1402 for (i = 0; i < nz; ++i)
1403 mat->cval[i] -= beg_row;
1407 #define __FUNC__ "Factor_dhMaxPivotInverse"
1409 Factor_dhMaxPivotInverse (
Factor_dh mat)
1411 START_FUNC_DH
int i, m = mat->m, *diags = mat->diag;
1412 REAL_DH *aval = mat->aval;
1413 double minGlobal = 0.0, min = aval[diags[0]];
1416 for (i = 0; i < m; ++i)
1417 min = MIN (min, fabs (aval[diags[i]]));
1424 MPI_Reduce (&min, &minGlobal, 1, MPI_DOUBLE, MPI_MIN, 0, comm_dh);
1433 retval = 1.0 / minGlobal;
1435 END_FUNC_VAL (retval)}
1438 #define __FUNC__ "Factor_dhMaxValue"
1442 START_FUNC_DH
double maxGlobal = 0.0, max = 0.0;
1443 int i, nz = mat->rp[mat->m];
1444 REAL_DH *aval = mat->aval;
1446 for (i = 0; i < nz; ++i)
1448 max = MAX (max, fabs (aval[i]));
1457 MPI_Reduce (&max, &maxGlobal, 1, MPI_DOUBLE, MPI_MAX, 0, comm_dh);
1459 END_FUNC_VAL (maxGlobal)}
1463 #define __FUNC__ "Factor_dhCondEst"
1467 START_FUNC_DH
double max = 0.0, maxGlobal = 0.0;
1472 Vec_dhCreate (&lhs);
1474 Vec_dhInit (lhs, m);
1476 Vec_dhDuplicate (lhs, &rhs);
1478 Vec_dhSet (rhs, 1.0);
1480 Euclid_dhApply (ctx, rhs->vals, lhs->vals);
1484 for (i = 0; i < m; ++i)
1486 max = MAX (max, fabs (x[i]));
1495 MPI_Reduce (&max, &maxGlobal, 1, MPI_DOUBLE, MPI_MAX, 0, comm_dh);
1497 END_FUNC_VAL (maxGlobal)}