68 #define __FUNC__ "Factor_dhCreate" 76 SET_V_ERROR (
"you must change MAX_MPI_TASKS and recompile!");
114 #define __FUNC__ "Factor_dhDestroy" 186 #define __FUNC__ "create_fake_mat_private" 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" 217 #define __FUNC__ "Factor_dhReadNz" 223 ierr = MPI_Allreduce (&nz, &retval, 1, MPI_INT, MPI_SUM,
comm_dh);
230 #define __FUNC__ "Factor_dhPrintRows" 235 int m = mat->
m, i, j;
249 "\n----------------------- Factor_dhPrintRows ------------------\n");
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]);
282 #define __FUNC__ "Factor_dhPrintDiags" 292 "\n----------------------- Factor_dhPrintDiags ------------------\n");
295 for (pe = 0; pe <
np_dh; ++pe)
300 fprintf (fp,
"----- subdomain: %i processor: %i\n", pe,
myid_dh);
301 for (i = 0; i <
m; ++i)
310 fprintf (fp,
"%i %g ZERO\n", i + 1 +
beg_row,
320 #define __FUNC__ "Factor_dhPrintGraph" 328 SET_V_ERROR (
"only implemented for single mpi task");
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" 370 int m = mat->
m, *
rp = mat->
rp;
386 for (pe = 0; pe <
np_dh; ++pe)
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,
446 #define __FUNC__ "setup_receives_private" 449 double *recvBuf, MPI_Request * req,
450 int *reqind,
int reqlen,
int *outlist,
bool debug)
458 "\nFACT ========================================================\n");
459 fprintf (
logFile,
"FACT STARTING: setup_receives_private\n");
462 for (i = 0; i < reqlen; i = j)
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]);
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,
514 #define __FUNC__ "setup_sends_private" 517 int *o2n_subdomain,
bool debug)
521 MPI_Status *statuses = mat->
status;
525 int myidNEW = o2n_subdomain[
myid_dh];
530 fprintf (
logFile,
"FACT \nSTARTING: setup_sends_private\n");
535 for (i = 0; i <
np_dh; i++)
539 if (o2n_subdomain[i] < myidNEW)
565 for (i = 0; i <
np_dh; i++)
569 isHigher = (o2n_subdomain[i] < myidNEW) ?
false :
true;
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;
628 fprintf (
logFile,
"FACT send to P_%i: ", i);
629 for (j = 0; j < inlist[i]; ++j)
630 fprintf (
logFile,
"%i ", rcvBuf[j] + 1);
647 #define __FUNC__ "Factor_dhSolveSetup" 670 for (i = 0; i <
np_dh; ++i)
674 end_rows[i] = beg_rows[i] + row_count[i];
690 fprintf (stderr,
"Numbering_dhSetup completed\n");
701 fprintf (
logFile,
"FACT num_extLo= %i num_extHi= %i\n",
730 MPI_Alltoall (outlist, 1, MPI_INT, inlist, 1, MPI_INT,
comm_dh);
741 for (row = 0; row <
m; row++)
743 int len =
rp[row + 1] -
rp[row];
744 int *ind =
cval +
rp[row];
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)
782 double *
aval,
double *rhs,
double *work_y,
787 double *
aval,
double *work_y,
788 double *work_x,
bool debug);
794 #define __FUNC__ "Factor_dhSolve" 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;
808 double *work_y = mat->work_y_lo;
809 double *work_x = mat->work_x_hi;
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);
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]);
896 MPI_Startall (mat->num_sendHi, mat->send_reqHi);
902 "\nFACT sending 'y' values to higher neighbor:\nFACT ");
917 ierr = MPI_Waitall (mat->num_recvHi, mat->recv_reqHi, mat->status);
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]);
938 work_y, work_x,
debug);
953 ierr = MPI_Startall (mat->num_sendLo, mat->send_reqLo);
960 "\nFACT sending 'x' values to lower neighbor:\nFACT ");
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]);
995 ierr = MPI_Waitall (mat->num_sendLo, mat->send_reqLo, mat->status);
1001 ierr = MPI_Waitall (mat->num_sendHi, mat->send_reqHi, mat->status);
1009 #define __FUNC__ "forward_solve_private" 1013 double *rhs,
double *work_y,
bool debug)
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];
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]);
1058 fprintf (
logFile,
"-----------\n");
1061 fprintf (
logFile,
"\nFACT work vector at end of forward solve:\n");
1062 for (i = 0; i < to; i++)
1068 for (i = from; i < to; ++i)
1070 int len =
diag[i] -
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" 1090 double *work_y,
double *work_x,
bool debug)
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;
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]);
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;
1129 double sum = work_y[i];
1131 for (j = 0; j < len; ++j)
1134 sum -= (val[j] * work_x[idx]);
1142 #define __FUNC__ "Factor_dhInit" 1145 double rho,
int id,
int beg_rowP,
Factor_dh * Fout)
1183 #define __FUNC__ "Factor_dhReallocate" 1189 if (used + additional > F->
alloc)
1192 while (
alloc < used + additional)
1198 memcpy (F->
cval, tmpI, used * sizeof (
int));
1206 memcpy (F->
fill, tmpI, used * sizeof (
int));
1223 #define __FUNC__ "Factor_dhTranspose" 1238 if (
B->aval ==
NULL)
1255 #define __FUNC__ "Factor_dhSolveSeq" 1262 printf (
"F is null.\n");
1266 printf (
"F isn't null.\n");
1269 int i, j, *vi, nz,
m = F->m;
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++)
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--)
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]);
1353 for (i = 1; i <
m; i++)
1360 sum -= (*v++ * work[*vi++]);
1365 for (i =
m - 1; i >= 0; i--)
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" 1388 int nz = mat->
rp[mat->
m];
1390 for (i = 0; i < nz; ++i)
1395 #define __FUNC__ "unadjust_bj_private" 1400 int nz = mat->
rp[mat->
m];
1402 for (i = 0; i < nz; ++i)
1407 #define __FUNC__ "Factor_dhMaxPivotInverse" 1413 double minGlobal = 0.0,
min =
aval[diags[0]];
1416 for (i = 0; i <
m; ++i)
1424 MPI_Reduce (&
min, &minGlobal, 1, MPI_DOUBLE, MPI_MIN, 0,
comm_dh);
1433 retval = 1.0 / minGlobal;
1438 #define __FUNC__ "Factor_dhMaxValue" 1443 int i, nz = mat->
rp[mat->
m];
1446 for (i = 0; i < nz; ++i)
1457 MPI_Reduce (&
max, &maxGlobal, 1, MPI_DOUBLE, MPI_MAX, 0,
comm_dh);
1463 #define __FUNC__ "Factor_dhCondEst" 1484 for (i = 0; i <
m; ++i)
1495 MPI_Reduce (&
max, &maxGlobal, 1, MPI_DOUBLE, MPI_MAX, 0,
comm_dh);
void Factor_dhPrintRows(Factor_dh mat, FILE *fp)
MPI_Request send_reqHi[MAX_MPI_TASKS]
static void setup_sends_private(Factor_dh mat, int *inlist, int *o2n_subdomain, bool debug)
MPI_Request recv_reqHi[MAX_MPI_TASKS]
static void destroy_fake_mat_private(Mat_dh matFake)
void Vec_dhSet(Vec_dh v, double value)
void Factor_dhPrintGraph(Factor_dh mat, char *filename)
void Vec_dhDuplicate(Vec_dh v, Vec_dh *out)
double Factor_dhCondEst(Factor_dh mat, Euclid_dh ctx)
void Factor_dhCreate(Factor_dh *mat)
#define END_FUNC_VAL(retval)
void Factor_dhInit(void *A, bool fillFlag, bool avalFlag, double rho, int id, int beg_rowP, Factor_dh *Fout)
void Euclid_dhApply(Euclid_dh ctx, double *rhs, double *lhs)
void Vec_dhInit(Vec_dh v, int size)
bool Parser_dhHasSwitch(Parser_dh p, char *s)
void Numbering_dhSetup(Numbering_dh numb, Mat_dh mat)
static void unadjust_bj_private(Factor_dh mat)
void Factor_dhTranspose(Factor_dh A, Factor_dh *Bout)
int mat_find_owner(int *beg_rows, int *end_rows, int index)
void Numbering_dhGlobalToLocal(Numbering_dh numb, int len, int *global, int *local)
void Factor_dhSolve(double *rhs, double *lhs, Euclid_dh ctx)
void Factor_dhPrintDiags(Factor_dh mat, FILE *fp)
#define CHECK_MPI_ERROR(errCode)
void fprintf_dh(FILE *fp, char *fmt,...)
void Factor_dhSolveSetup(Factor_dh mat, SubdomainGraph_dh sg)
void Mat_dhCreate(Mat_dh *mat)
double Factor_dhMaxValue(Factor_dh mat)
void mat_dh_transpose_private(int m, int *RP, int **rpOUT, int *CVAL, int **cvalOUT, double *AVAL, double **avalOUT)
MPI_Request send_reqLo[MAX_MPI_TASKS]
#define CHECK_ERROR(retval)
void Numbering_dhCreate(Numbering_dh *numb)
MPI_Status status[MAX_MPI_TASKS]
void closeFile_dh(FILE *fpIN)
void Numbering_dhDestroy(Numbering_dh numb)
static void forward_solve_private(int m, int from, int to, int *rp, int *cval, int *diag, double *aval, double *rhs, double *work_y, bool debug)
void Vec_dhCreate(Vec_dh *v)
void Factor_dhPrintTriples(Factor_dh mat, char *filename)
static void backward_solve_private(int m, int from, int to, int *rp, int *cval, int *diag, double *aval, double *work_y, double *work_x, bool debug)
MPI_Request recv_reqLo[MAX_MPI_TASKS]
static int setup_receives_private(Factor_dh mat, int *beg_rows, int *end_rows, double *recvBuf, MPI_Request *req, int *reqind, int reqlen, int *outlist, bool debug)
void Factor_dhReallocate(Factor_dh F, int used, int additional)
MPI_Request requests[MAX_MPI_TASKS]
static void create_fake_mat_private(Factor_dh mat, Mat_dh *matFakeIN)
FILE * openFile_dh(const char *filenameIN, const char *modeIN)
void Mat_dhDestroy(Mat_dh mat)
double Factor_dhMaxPivotInverse(Factor_dh mat)
void EuclidGetDimensions(void *A, int *beg_row, int *rowsLocal, int *rowsGlobal)
int Factor_dhReadNz(Factor_dh mat)
#define CHECK_MPI_V_ERROR(errCode)
void Factor_dhSolveSeq(double *rhs, double *lhs, Euclid_dh ctx)
void Factor_dhDestroy(Factor_dh mat)
static void adjust_bj_private(Factor_dh mat)