44 #include "getRow_dh.h"
45 #include "SubdomainGraph_dh.h"
46 #include "TimeLog_dh.h"
48 #include "Numbering_dh.h"
49 #include "Parser_dh.h"
50 #include "mat_dh_private.h"
52 #include "Hash_i_dh.h"
54 static void setup_matvec_sends_private (
Mat_dh mat,
int *inlist);
55 static void setup_matvec_receives_private (
Mat_dh mat,
int *beg_rows,
56 int *end_rows,
int reqlen,
57 int *reqind,
int *outlist);
63 not used anyplace, I think;
65 expansion ?[mar 21, 2 K + 1]
66 static
void Mat_dhAllocate_getRow_private (
Mat_dh A);
69 static bool commsOnly =
false;
72 #define __FUNC__ "Mat_dhCreate"
73 void Mat_dhCreate (
Mat_dh * mat)
81 commsOnly = Parser_dhHasSwitch (parser_dh,
"-commsOnly");
82 if (myid_dh == 0 && commsOnly ==
true)
101 tmp->len_private = 0;
102 tmp->rowCheckedOut = -1;
103 tmp->cval_private = NULL;
104 tmp->aval_private = NULL;
106 tmp->row_perm = NULL;
110 tmp->recv_req = NULL;
111 tmp->send_req = NULL;
119 tmp->matvecIsSetup =
false;
121 Mat_dhZeroTiming (tmp);
123 tmp->matvec_timing =
true;
125 tmp->debug = Parser_dhHasSwitch (parser_dh,
"-debug_Mat");
129 #define __FUNC__ "Mat_dhDestroy"
131 Mat_dhDestroy (
Mat_dh mat)
142 if (mat->len != NULL)
147 if (mat->cval != NULL)
152 if (mat->aval != NULL)
157 if (mat->diag != NULL)
162 if (mat->fill != NULL)
167 if (mat->cval_private != NULL)
169 FREE_DH (mat->cval_private);
172 if (mat->aval_private != NULL)
174 FREE_DH (mat->aval_private);
177 if (mat->row_perm != NULL)
179 FREE_DH (mat->row_perm);
184 for (i = 0; i < mat->num_recv; i++)
185 MPI_Request_free (&mat->recv_req[i]);
186 for (i = 0; i < mat->num_send; i++)
187 MPI_Request_free (&mat->send_req[i]);
188 if (mat->recv_req != NULL)
190 FREE_DH (mat->recv_req);
193 if (mat->send_req != NULL)
195 FREE_DH (mat->send_req);
198 if (mat->status != NULL)
200 FREE_DH (mat->status);
203 if (mat->recvbuf != NULL)
205 FREE_DH (mat->recvbuf);
208 if (mat->sendbuf != NULL)
210 FREE_DH (mat->sendbuf);
213 if (mat->sendind != NULL)
215 FREE_DH (mat->sendind);
219 if (mat->matvecIsSetup)
221 Mat_dhMatVecSetdown (mat);
224 if (mat->numb != NULL)
226 Numbering_dhDestroy (mat->numb);
236 #define __FUNC__ "Mat_dhMatVecSetDown"
238 Mat_dhMatVecSetdown (
Mat_dh mat)
240 START_FUNC_DH
if (ignoreMe)
241 SET_V_ERROR (
"not implemented");
247 #define __FUNC__ "Mat_dhMatVecSetup"
249 Mat_dhMatVecSetup (
Mat_dh mat)
251 START_FUNC_DH
if (np_dh == 1)
258 int *outlist, *inlist;
259 int ierr, i, row, *rp = mat->rp, *cval = mat->cval;
262 int firstLocal = mat->beg_row;
263 int lastLocal = firstLocal + m;
264 int *beg_rows, *end_rows;
267 (MPI_Request *) MALLOC_DH (np_dh *
sizeof (MPI_Request));
270 (MPI_Request *) MALLOC_DH (np_dh *
sizeof (MPI_Request));
272 mat->status = (MPI_Status *) MALLOC_DH (np_dh *
sizeof (MPI_Status));
274 beg_rows = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
276 end_rows = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
287 MPI_Allgather (&firstLocal, 1, MPI_INT, beg_rows, 1, MPI_INT,
290 CHECK_MPI_V_ERROR (ierr);
293 MPI_Allgather (&lastLocal, 1, MPI_INT, end_rows, 1, MPI_INT,
295 CHECK_MPI_V_ERROR (ierr);
298 outlist = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
300 inlist = (
int *) MALLOC_DH (np_dh *
sizeof (
int));
302 for (i = 0; i < np_dh; ++i)
309 Numbering_dhCreate (&(mat->numb));
312 Numbering_dhSetup (numb, mat);
315 setup_matvec_receives_private (mat, beg_rows, end_rows, numb->num_ext,
316 numb->idx_ext, outlist);
321 inlist[0] = outlist[0];
326 MPI_Alltoall (outlist, 1, MPI_INT, inlist, 1, MPI_INT, comm_dh);
327 CHECK_MPI_V_ERROR (ierr);
330 setup_matvec_sends_private (mat, inlist);
334 for (row = 0; row < m; row++)
336 int len = rp[row + 1] - rp[row];
337 int *ind = cval + rp[row];
338 Numbering_dhGlobalToLocal (numb, len, ind, ind);
358 #define __FUNC__ "setup_matvec_receives_private"
360 setup_matvec_receives_private (
Mat_dh mat,
int *beg_rows,
int *end_rows,
361 int reqlen,
int *reqind,
int *outlist)
363 START_FUNC_DH
int ierr, i, j, this_pe;
371 mat->recvbuf = (
double *) MALLOC_DH ((reqlen + m) *
sizeof (double));
373 for (i = 0; i < reqlen; i = j)
376 this_pe = mat_find_owner (beg_rows, end_rows, reqind[i]);
380 for (j = i + 1; j < reqlen; j++)
383 if (reqind[j] < beg_rows[this_pe] || reqind[j] > end_rows[this_pe])
389 MPI_Isend (&reqind[i], j - i, MPI_INT, this_pe, 444, comm_dh,
391 CHECK_MPI_V_ERROR (ierr);
392 ierr = MPI_Request_free (&request);
393 CHECK_MPI_V_ERROR (ierr);
396 outlist[this_pe] = j - i;
399 MPI_Recv_init (&mat->recvbuf[i + m], j - i, MPI_DOUBLE, this_pe, 555,
400 comm_dh, &mat->recv_req[mat->num_recv]);
401 CHECK_MPI_V_ERROR (ierr);
404 mat->recvlen += j - i;
411 #define __FUNC__ "setup_matvec_sends_private"
413 setup_matvec_sends_private (
Mat_dh mat,
int *inlist)
415 START_FUNC_DH
int ierr, i, j, sendlen, first = mat->beg_row;
416 MPI_Request *requests;
417 MPI_Status *statuses;
419 requests = (MPI_Request *) MALLOC_DH (np_dh *
sizeof (MPI_Request));
421 statuses = (MPI_Status *) MALLOC_DH (np_dh *
sizeof (MPI_Status));
426 for (i = 0; i < np_dh; i++)
427 sendlen += inlist[i];
428 mat->sendlen = sendlen;
429 mat->sendbuf = (
double *) MALLOC_DH (sendlen *
sizeof (
double));
431 mat->sendind = (
int *) MALLOC_DH (sendlen *
sizeof (
int));
436 for (i = 0; i < np_dh; i++)
442 MPI_Irecv (&mat->sendind[j], inlist[i], MPI_INT, i, 444, comm_dh,
443 &requests[mat->num_send]);
444 CHECK_MPI_V_ERROR (ierr);
447 MPI_Send_init (&mat->sendbuf[j], inlist[i], MPI_DOUBLE, i, 555,
448 comm_dh, &mat->send_req[mat->num_send]);
449 CHECK_MPI_V_ERROR (ierr);
457 mat->time[MATVEC_WORDS] = j;
460 ierr = MPI_Waitall (mat->num_send, requests, statuses);
461 CHECK_MPI_V_ERROR (ierr);
464 for (i = 0; i < mat->sendlen; i++)
465 mat->sendind[i] -= first;
474 #define __FUNC__ "Mat_dhMatVec"
476 Mat_dhMatVec (
Mat_dh mat,
double *x,
double *b)
478 START_FUNC_DH
if (np_dh == 1)
480 Mat_dhMatVec_uni (mat, x, b);
486 int ierr, i, row, m = mat->m;
487 int *rp = mat->rp, *cval = mat->cval;
488 double *aval = mat->aval;
489 int *sendind = mat->sendind;
490 int sendlen = mat->sendlen;
491 double *sendbuf = mat->sendbuf;
492 double *recvbuf = mat->recvbuf;
493 double t1 = 0, t2 = 0, t3 = 0, t4 = 0;
494 bool timeFlag = mat->matvec_timing;
503 for (i = 0; i < sendlen; i++)
504 sendbuf[i] = x[sendind[i]];
510 mat->time[MATVEC_TIME] += (t2 - t1);
514 ierr = MPI_Startall (mat->num_recv, mat->recv_req);
515 CHECK_MPI_V_ERROR (ierr);
516 ierr = MPI_Startall (mat->num_send, mat->send_req);
517 CHECK_MPI_V_ERROR (ierr);
518 ierr = MPI_Waitall (mat->num_recv, mat->recv_req, mat->status);
519 CHECK_MPI_V_ERROR (ierr);
520 ierr = MPI_Waitall (mat->num_send, mat->send_req, mat->status);
521 CHECK_MPI_V_ERROR (ierr);
527 mat->time[MATVEC_MPI_TIME] += (t3 - t2);
533 for (i = 0; i < m; i++)
537 for (row = 0; row < m; row++)
539 int len = rp[row + 1] - rp[row];
540 int *ind = cval + rp[row];
541 double *val = aval + rp[row];
543 for (i = 0; i < len; i++)
545 temp += (val[i] * recvbuf[ind[i]]);
554 mat->time[MATVEC_TOTAL_TIME] += (t4 - t1);
555 mat->time[MATVEC_TIME] += (t4 - t3);
562 #define __FUNC__ "Mat_dhMatVec_omp"
564 Mat_dhMatVec_omp (
Mat_dh mat,
double *x,
double *b)
566 START_FUNC_DH
int ierr, i, row, m = mat->m;
567 int *rp = mat->rp, *cval = mat->cval;
568 double *aval = mat->aval;
569 int *sendind = mat->sendind;
570 int sendlen = mat->sendlen;
571 double *sendbuf = mat->sendbuf;
572 double *recvbuf = mat->recvbuf;
573 double t1 = 0, t2 = 0, t3 = 0, t4 = 0, tx = 0;
576 bool timeFlag = mat->matvec_timing;
582 for (i = 0; i < sendlen; i++)
583 sendbuf[i] = x[sendind[i]];
588 mat->time[MATVEC_TIME] += (t2 - t1);
591 ierr = MPI_Startall (mat->num_recv, mat->recv_req);
592 CHECK_MPI_V_ERROR (ierr);
593 ierr = MPI_Startall (mat->num_send, mat->send_req);
594 CHECK_MPI_V_ERROR (ierr);
595 ierr = MPI_Waitall (mat->num_recv, mat->recv_req, mat->status);
596 CHECK_MPI_V_ERROR (ierr);
597 ierr = MPI_Waitall (mat->num_send, mat->send_req, mat->status);
598 CHECK_MPI_V_ERROR (ierr);
603 mat->time[MATVEC_MPI_TIME] += (t3 - t2);
607 for (i = 0; i < m; i++)
613 mat->time[MATVEC_MPI_TIME2] += (tx - t1);
618 for (row = 0; row < m; row++)
620 len = rp[row + 1] - rp[row];
621 ind = cval + rp[row];
622 val = aval + rp[row];
624 for (i = 0; i < len; i++)
626 temp += (val[i] * recvbuf[ind[i]]);
634 mat->time[MATVEC_TOTAL_TIME] += (t4 - t1);
635 mat->time[MATVEC_TIME] += (t4 - t3);
643 #define __FUNC__ "Mat_dhMatVec_uni_omp"
645 Mat_dhMatVec_uni_omp (
Mat_dh mat,
double *x,
double *b)
647 START_FUNC_DH
int i, row, m = mat->m;
648 int *rp = mat->rp, *cval = mat->cval;
649 double *aval = mat->aval;
650 double t1 = 0, t2 = 0;
651 bool timeFlag = mat->matvec_timing;
659 for (row = 0; row < m; row++)
661 int len = rp[row + 1] - rp[row];
662 int *ind = cval + rp[row];
663 double *val = aval + rp[row];
665 for (i = 0; i < len; i++)
667 temp += (val[i] * x[ind[i]]);
675 mat->time[MATVEC_TIME] += (t2 - t1);
676 mat->time[MATVEC_TOTAL_TIME] += (t2 - t1);
684 #define __FUNC__ "Mat_dhMatVec_uni"
686 Mat_dhMatVec_uni (
Mat_dh mat,
double *x,
double *b)
688 START_FUNC_DH
int i, row, m = mat->m;
689 int *rp = mat->rp, *cval = mat->cval;
690 double *aval = mat->aval;
691 double t1 = 0, t2 = 0;
692 bool timeFlag = mat->matvec_timing;
697 for (row = 0; row < m; row++)
699 int len = rp[row + 1] - rp[row];
700 int *ind = cval + rp[row];
701 double *val = aval + rp[row];
703 for (i = 0; i < len; i++)
705 temp += (val[i] * x[ind[i]]);
713 mat->time[MATVEC_TIME] += (t2 - t1);
714 mat->time[MATVEC_TOTAL_TIME] += (t2 - t1);
721 #define __FUNC__ "Mat_dhReadNz"
725 START_FUNC_DH
int ierr, retval = mat->rp[mat->m];
727 ierr = MPI_Allreduce (&nz, &retval, 1, MPI_INT, MPI_SUM, comm_dh);
728 CHECK_MPI_ERROR (ierr);
729 END_FUNC_VAL (retval)}
736 #define __FUNC__ "Mat_dhAllocate_getRow_private"
738 Mat_dhAllocate_getRow_private (
Mat_dh A)
740 START_FUNC_DH
int i, *rp = A->rp, len = 0;
744 for (i = 0; i < m; ++i)
745 len = MAX (len, rp[i + 1] - rp[i]);
749 if (len > A->len_private)
751 if (A->cval_private != NULL)
753 FREE_DH (A->cval_private);
756 if (A->aval_private != NULL)
758 FREE_DH (A->aval_private);
764 A->cval_private = (
int *) MALLOC_DH (len *
sizeof (
int));
766 A->aval_private = (
double *) MALLOC_DH (len *
sizeof (
double));
768 A->len_private = len;
774 #define __FUNC__ "Mat_dhZeroTiming"
776 Mat_dhZeroTiming (
Mat_dh mat)
780 for (i = 0; i < MAT_DH_BINS; ++i)
783 mat->time_max[i] = 0;
784 mat->time_min[i] = 0;
789 #define __FUNC__ "Mat_dhReduceTiming"
791 Mat_dhReduceTiming (
Mat_dh mat)
793 START_FUNC_DH
if (mat->time[MATVEC_MPI_TIME])
795 mat->time[MATVEC_RATIO] =
796 mat->time[MATVEC_TIME] / mat->time[MATVEC_MPI_TIME];
798 MPI_Allreduce (mat->time, mat->time_min, MAT_DH_BINS, MPI_DOUBLE, MPI_MIN,
800 MPI_Allreduce (mat->time, mat->time_max, MAT_DH_BINS, MPI_DOUBLE, MPI_MAX,
805 #define __FUNC__ "Mat_dhPermute"
810 int i, j, *RP = A->rp, *CVAL = A->cval;
811 int *o2n, *rp, *cval, m = A->m, nz = RP[m];
812 double *aval, *AVAL = A->aval;
820 o2n = (
int *) MALLOC_DH (m *
sizeof (
int));
822 for (i = 0; i < m; ++i)
826 rp = B->rp = (
int *) MALLOC_DH ((m + 1) *
sizeof (int));
828 cval = B->cval = (
int *) MALLOC_DH (nz *
sizeof (
int));
830 aval = B->aval = (
double *) MALLOC_DH (nz *
sizeof (
double));
835 for (i = 0; i < m; ++i)
838 rp[i + 1] = RP[oldRow + 1] - RP[oldRow];
840 for (i = 1; i <= m; ++i)
841 rp[i] = rp[i] + rp[i - 1];
843 for (i = 0; i < m; ++i)
847 for (j = RP[oldRow]; j < RP[oldRow + 1]; ++j)
849 cval[idx] = o2n[CVAL[j]];
866 #define __FUNC__ "Mat_dhPrintGraph"
870 START_FUNC_DH
int pe,
id = myid_dh;
875 id = sg->o2n_sub[id];
878 for (pe = 0; pe < np_dh; ++pe)
880 ierr = MPI_Barrier (comm_dh);
881 CHECK_MPI_V_ERROR (ierr);
886 mat_dh_print_graph_private (A->m, A->beg_row, A->rp, A->cval,
887 A->aval, NULL, NULL, NULL, fp);
892 int beg_row = sg->beg_rowP[myid_dh];
893 mat_dh_print_graph_private (A->m, beg_row, A->rp, A->cval,
894 A->aval, sg->n2o_row, sg->o2n_col,
904 #define __FUNC__ "Mat_dhPrintRows"
908 START_FUNC_DH
bool noValues;
909 int m = A->m, *rp = A->rp, *cval = A->cval;
910 double *aval = A->aval;
912 noValues = (Parser_dhHasSwitch (parser_dh,
"-noValues"));
922 int beg_row = A->beg_row;
925 "\n----- A, unpermuted ------------------------------------\n");
926 for (i = 0; i < m; ++i)
928 fprintf (fp,
"%i :: ", 1 + i + beg_row);
929 for (j = rp[i]; j < rp[i + 1]; ++j)
933 fprintf (fp,
"%i ", 1 + cval[j]);
937 fprintf (fp,
"%i,%g ; ", 1 + cval[j], aval[j]);
952 for (i = 0; i < sg->blocks; ++i)
954 int oldBlock = sg->n2o_sub[i];
959 int beg_row = sg->beg_row[oldBlock];
960 int end_row = beg_row + sg->row_count[oldBlock];
964 "\n----- A, permuted, single mpi task ------------------\n");
965 fprintf (fp,
"---- new subdomain: %i; old subdomain: %i\n", i,
967 fprintf (fp,
" old beg_row: %i; new beg_row: %i\n",
968 sg->beg_row[oldBlock], sg->beg_rowP[oldBlock]);
969 fprintf (fp,
" local rows in this block: %i\n",
970 sg->row_count[oldBlock]);
971 fprintf (fp,
" bdry rows in this block: %i\n",
972 sg->bdry_count[oldBlock]);
973 fprintf (fp,
" 1st bdry row= %i \n",
974 1 + end_row - sg->bdry_count[oldBlock]);
976 for (oldRow = beg_row; oldRow < end_row; ++oldRow)
981 fprintf (fp,
"%3i (old= %3i) :: ", idx, 1 + oldRow);
983 Mat_dhGetRow (A, oldRow, &len, &cval, &aval);
986 for (k = 0; k < len; ++k)
990 fprintf (fp,
"%i ", 1 + sg->o2n_col[cval[k]]);
994 fprintf (fp,
"%i,%g ; ", 1 + sg->o2n_col[cval[k]],
1000 Mat_dhRestoreRow (A, oldRow, &len, &cval, &aval);
1012 int *o2n_col = sg->o2n_col, *n2o_row = sg->n2o_row;
1013 int beg_row = sg->beg_row[myid_dh];
1014 int beg_rowP = sg->beg_rowP[myid_dh];
1017 for (i = 0; i < m; ++i)
1019 int row = n2o_row[i];
1020 fprintf (fp,
"%3i (old= %3i) :: ", 1 + i + beg_rowP,
1022 for (j = rp[row]; j < rp[row + 1]; ++j)
1028 if (col >= beg_row && col < beg_row + m)
1030 col = o2n_col[col - beg_row] + beg_rowP;
1037 tmp = Hash_i_dhLookup (hash, col);
1042 "nonlocal column= %i not in hash table",
1044 SET_V_ERROR (msgBuf_dh);
1054 fprintf (fp,
"%i ", 1 + col);
1058 fprintf (fp,
"%i,%g ; ", 1 + col, aval[j]);
1069 #define __FUNC__ "Mat_dhPrintTriples"
1073 START_FUNC_DH
int m = A->m, *rp = A->rp, *cval = A->cval;
1074 double *aval = A->aval;
1079 noValues = (Parser_dhHasSwitch (parser_dh,
"-noValues"));
1082 matlab = (Parser_dhHasSwitch (parser_dh,
"-matlab"));
1090 int beg_row = A->beg_row;
1093 for (pe = 0; pe < np_dh; ++pe)
1095 MPI_Barrier (comm_dh);
1100 fp = openFile_dh (filename,
"w");
1105 fp = openFile_dh (filename,
"a");
1109 for (i = 0; i < m; ++i)
1111 for (j = rp[i]; j < rp[i + 1]; ++j)
1115 fprintf (fp,
"%i %i\n", 1 + i + beg_row,
1121 if (val == 0.0 && matlab)
1122 val = _MATLAB_ZERO_;
1123 fprintf (fp, TRIPLES_FORMAT, 1 + i + beg_row,
1137 else if (np_dh == 1)
1139 int i, j, k, idx = 1;
1141 fp = openFile_dh (filename,
"w");
1144 for (i = 0; i < sg->blocks; ++i)
1146 int oldBlock = sg->n2o_sub[i];
1147 int beg_row = sg->beg_rowP[oldBlock];
1148 int end_row = beg_row + sg->row_count[oldBlock];
1150 for (j = beg_row; j < end_row; ++j)
1154 int oldRow = sg->n2o_row[j];
1156 Mat_dhGetRow (A, oldRow, &len, &cval, &aval);
1161 for (k = 0; k < len; ++k)
1163 fprintf (fp,
"%i %i\n", idx, 1 + sg->o2n_col[cval[k]]);
1170 for (k = 0; k < len; ++k)
1172 double val = aval[k];
1173 if (val == 0.0 && matlab)
1174 val = _MATLAB_ZERO_;
1175 fprintf (fp, TRIPLES_FORMAT, idx,
1176 1 + sg->o2n_col[cval[k]], val);
1180 Mat_dhRestoreRow (A, oldRow, &len, &cval, &aval);
1192 int *o2n_col = sg->o2n_col, *n2o_row = sg->n2o_row;
1193 int beg_row = sg->beg_row[myid_dh];
1194 int beg_rowP = sg->beg_rowP[myid_dh];
1196 int id = sg->o2n_sub[myid_dh];
1198 for (pe = 0; pe < np_dh; ++pe)
1200 MPI_Barrier (comm_dh);
1205 fp = openFile_dh (filename,
"w");
1210 fp = openFile_dh (filename,
"a");
1214 for (i = 0; i < m; ++i)
1216 int row = n2o_row[i];
1217 for (j = rp[row]; j < rp[row + 1]; ++j)
1224 if (val == 0.0 && matlab)
1225 val = _MATLAB_ZERO_;
1229 if (col >= beg_row && col < beg_row + m)
1231 col = o2n_col[col - beg_row] + beg_rowP;
1238 tmp = Hash_i_dhLookup (hash, col);
1243 "nonlocal column= %i not in hash table",
1245 SET_V_ERROR (msgBuf_dh);
1255 fprintf (fp,
"%i %i\n", 1 + i + beg_rowP, 1 + col);
1259 fprintf (fp, TRIPLES_FORMAT, 1 + i + beg_rowP,
1274 #define __FUNC__ "Mat_dhPrintCSR"
1278 START_FUNC_DH FILE * fp;
1282 SET_V_ERROR (
"only implemented for a single mpi task");
1287 (
"not implemented for reordered matrix (SubdomainGraph_dh should be NULL)");
1290 fp = openFile_dh (filename,
"w");
1295 mat_dh_print_csr_private (A->m, A->rp, A->cval, A->aval, fp);
1300 mat_dh_print_csr_private (A->m, A->rp, A->cval, A->aval, fp);
1310 #define __FUNC__ "Mat_dhPrintBIN"
1314 START_FUNC_DH
if (np_dh > 1)
1316 SET_V_ERROR (
"only implemented for a single MPI task");
1322 SET_V_ERROR (
"not implemented for reordering; ensure sg=NULL");
1325 io_dh_print_ebin_mat_private (A->m, A->beg_row, A->rp, A->cval, A->aval,
1326 NULL, NULL, NULL, filename);
1336 #define __FUNC__ "Mat_dhReadCSR"
1338 Mat_dhReadCSR (
Mat_dh * mat,
char *filename)
1345 SET_V_ERROR (
"only implemented for a single MPI task");
1348 fp = openFile_dh (filename,
"r");
1353 mat_dh_read_csr_private (&A->m, &A->rp, &A->cval, &A->aval, fp);
1364 #define __FUNC__ "Mat_dhReadTriples"
1366 Mat_dhReadTriples (
Mat_dh * mat,
int ignore,
char *filename)
1368 START_FUNC_DH FILE * fp = NULL;
1373 SET_V_ERROR (
"only implemented for a single MPI task");
1376 fp = openFile_dh (filename,
"r");
1381 mat_dh_read_triples_private (ignore, &A->m, &A->rp, &A->cval, &A->aval, fp);
1395 #define __FUNC__ "Mat_dhReadBIN"
1397 Mat_dhReadBIN (
Mat_dh * mat,
char *filename)
1403 SET_V_ERROR (
"only implemented for a single MPI task");
1408 io_dh_read_ebin_mat_private (&A->m, &A->rp, &A->cval, &A->aval, filename);
1415 #define __FUNC__ "Mat_dhTranspose"
1423 SET_V_ERROR (
"only for sequential");
1430 mat_dh_transpose_private (A->m, A->rp, &B->rp, A->cval, &B->cval,
1436 #define __FUNC__ "Mat_dhMakeStructurallySymmetric"
1438 Mat_dhMakeStructurallySymmetric (
Mat_dh A)
1440 START_FUNC_DH
if (np_dh > 1)
1442 SET_V_ERROR (
"only for sequential");
1444 make_symmetric_private (A->m, &A->rp, &A->cval, &A->aval);
1448 void insert_diags_private (
Mat_dh A,
int ct);
1455 #define __FUNC__ "Mat_dhFixDiags"
1457 Mat_dhFixDiags (
Mat_dh A)
1459 START_FUNC_DH
int i, j;
1460 int *rp = A->rp, *cval = A->cval, m = A->m;
1462 double *aval = A->aval;
1465 for (i = 0; i < m; ++i)
1468 for (j = rp[i]; j < rp[i + 1]; ++j)
1485 (
"\nMat_dhFixDiags:: %i diags not explicitly present; inserting!\n",
1487 insert_diags_private (A, ct);
1495 for (i = 0; i < m; ++i)
1498 for (j = rp[i]; j < rp[i + 1]; ++j)
1500 sum += fabs (aval[j]);
1502 for (j = rp[i]; j < rp[i + 1]; ++j)
1514 #define __FUNC__ "insert_diags_private"
1516 insert_diags_private (
Mat_dh A,
int ct)
1518 START_FUNC_DH
int *RP = A->rp, *CVAL = A->cval;
1519 int *rp, *cval, m = A->m;
1520 double *aval, *AVAL = A->aval;
1521 int nz = RP[m] + ct;
1524 rp = A->rp = (
int *) MALLOC_DH ((m + 1) *
sizeof (int));
1526 cval = A->cval = (
int *) MALLOC_DH (nz *
sizeof (
int));
1528 aval = A->aval = (
double *) MALLOC_DH (nz *
sizeof (
double));
1532 for (i = 0; i < m; ++i)
1535 for (j = RP[i]; j < RP[i + 1]; ++j)
1537 cval[idx] = CVAL[j];
1538 aval[idx] = AVAL[j];
1563 #define __FUNC__ "Mat_dhPrintDiags"
1565 Mat_dhPrintDiags (
Mat_dh A, FILE * fp)
1567 START_FUNC_DH
int i, j, m = A->m;
1568 int *rp = A->rp, *cval = A->cval;
1569 double *aval = A->aval;
1572 "=================== diagonal elements ====================\n");
1573 for (i = 0; i < m; ++i)
1576 for (j = rp[i]; j < rp[i + 1]; ++j)
1580 fprintf (fp,
"%i %g\n", i + 1, aval[j]);
1587 fprintf (fp,
"%i ---------- missing\n", i + 1);
1594 #define __FUNC__ "Mat_dhGetRow"
1596 Mat_dhGetRow (
Mat_dh B,
int globalRow,
int *len,
int **ind,
double **val)
1598 START_FUNC_DH
int row = globalRow - B->beg_row;
1602 "requested globalRow= %i, which is local row= %i, but only have %i rows!",
1603 globalRow, row, B->m);
1604 SET_V_ERROR (msgBuf_dh);
1606 *len = B->rp[row + 1] - B->rp[row];
1608 *ind = B->cval + B->rp[row];
1610 *val = B->aval + B->rp[row];
1614 #define __FUNC__ "Mat_dhRestoreRow"
1616 Mat_dhRestoreRow (
Mat_dh B,
int row,
int *len,
int **ind,
double **val)
1618 START_FUNC_DH END_FUNC_DH}
1621 #define __FUNC__ "Mat_dhRowPermute"
1623 Mat_dhRowPermute (
Mat_dh mat)
1625 START_FUNC_DH
if (ignoreMe)
1626 SET_V_ERROR (
"turned off; compilation problem on blue");
1629 int i, j, m = mat->m, nz = mat->rp[m];
1633 bool debug = mat->debug;
1638 * = 1:Compute a row permutation of the matrix so that the
1639 * permuted matrix has as many entries on its diagonal as
1640 * possible.The values on the diagonal are of arbitrary size.
1641 * HSL subroutine MC21A / AD is used
for this
1642 . * = 2: Compute a row permutation of the matrix so that the smallest * value on the diagonal of the permuted matrix is maximized. * = 3:Compute a row permutation of the matrix so that the smallest
1643 * value on the diagonal of the permuted matrix is maximized.
1644 * The algorithm differs from the one used
for JOB
1645 = 2 and may * have quite a different performance. * = 4: Compute a row permutation of the matrix so that the sum * of the diagonal entries of the permuted matrix is maximized. * = 5:Compute a row permutation of the matrix so that the product
1646 * of the diagonal entries of the permuted matrix is maximized
1647 * and vectors to scale the matrix so that the nonzero diagonal
1648 * entries of the permuted matrix are one in absolute value and
1649 * all the off - diagonal entries are less than or equal to one in
1652 Parser_dhReadInt (parser_dh,
"-rowPermute", &algo);
1658 sprintf (msgBuf_dh,
"calling row permutation with algo= %i", algo);
1659 SET_INFO (msgBuf_dh);
1661 r1 = (
double *) MALLOC_DH (m *
sizeof (
double));
1663 c1 = (
double *) MALLOC_DH (m *
sizeof (
double));
1665 if (mat->row_perm == NULL)
1667 mat->row_perm = o2n = (
int *) MALLOC_DH (m *
sizeof (
int));
1672 o2n = mat->row_perm;
1675 Mat_dhTranspose (mat, &B);
1679 dldperm (algo, m, nz, B->rp, B->cval, B->aval, o2n, r1, c1);
1683 for (i = 0; i < nz; ++i)
1684 cval[i] = o2n[cval[i]];
1687 if (debug && logFile != NULL)
1689 fprintf (logFile,
"\n-------- row permutation vector --------\n");
1690 for (i = 0; i < m; ++i)
1691 fprintf (logFile,
"%i ", 1 + o2n[i]);
1692 fprintf (logFile,
"\n");
1696 printf (
"\n-------- row permutation vector --------\n");
1697 for (i = 0; i < m; ++i)
1698 printf (
"%i ", 1 + o2n[i]);
1705 for (i = 0; i < m; ++i)
1716 printf (
"@@@ [%i] Mat_dhRowPermute :: got natural ordering!\n",
1721 int *rp = B->rp, *cval = B->cval;
1722 double *aval = B->aval;
1727 (
"@@@ [%i] Mat_dhRowPermute :: scaling matrix rows and columns!\n",
1731 for (i = 0; i < m; i++)
1733 r1[i] = exp (r1[i]);
1734 c1[i] = exp (c1[i]);
1736 for (i = 0; i < m; i++)
1737 for (j = rp[i]; j < rp[i + 1]; j++)
1738 aval[j] *= r1[cval[j]] * c1[i];
1741 mat_dh_transpose_reuse_private (B->m, B->rp, B->cval, B->aval,
1742 mat->rp, mat->cval, mat->aval);
1760 #define __FUNC__ "Mat_dhPartition"
1762 build_adj_lists_private (
Mat_dh mat,
int **rpOUT,
int **cvalOUT)
1764 START_FUNC_DH
int m = mat->m;
1765 int *RP = mat->rp, *CVAL = mat->cval;
1767 int i, j, *rp, *cval, idx = 0;
1769 rp = *rpOUT = (
int *) MALLOC_DH ((m + 1) *
sizeof (int));
1771 cval = *cvalOUT = (
int *) MALLOC_DH (nz *
sizeof (
int));
1776 for (i = 0; i < m; ++i)
1778 for (j = RP[i]; j < RP[i + 1]; ++j)
1792 #define __FUNC__ "Mat_dhPartition"
1794 Mat_dhPartition (
Mat_dh mat,
int blocks,
1795 int **beg_rowOUT,
int **row_countOUT,
int **n2oOUT,
1799 #ifndef HAVE_METIS_DH
1801 SET_V_ERROR (
"not compiled for metis!");
1804 int *beg_row, *row_count, *n2o, *o2n, bk,
new, *part;
1806 int i, cutEdgeCount;
1808 int metisOpts[5] = { 0, 0, 0, 0, 0 };
1812 beg_row = *beg_rowOUT = (
int *) MALLOC_DH (blocks *
sizeof (
int));
1814 row_count = *row_countOUT = (
int *) MALLOC_DH (blocks *
sizeof (
int));
1816 *n2oOUT = n2o = (
int *) MALLOC_DH (m *
sizeof (
int));
1818 *o2nOUT = o2n = (
int *) MALLOC_DH (m *
sizeof (
int));
1822 == == == == == == == == == == == == == == == == == == == == == == == == == == == == == == = Metis arguments:
1824 n - number of nodes rp[], cval[]NULL, NULL, 0
1826 blocksIN, options[5] = 0::0 / 1 use defauls;
1831 == == == == == == == == == == == == == == == == == == == == == == == == ==
1835 build_adj_lists_private (mat, &rp, &cval);
1837 part = (
int *) MALLOC_DH (m *
sizeof (
int));
1841 METIS_PartGraphKway (&m, rp, cval, NULL, NULL,
1842 &zero, &zero, &blocks, metisOpts, &cutEdgeCount, part);
1851 printf_dh (
"\nmetis partitioning vector; blocks= %i\n", blocks);
1852 for (i = 0; i < m; ++i)
1853 printf_dh (
" %i %i\n", i + 1, part[i]);
1857 for (i = 0; i < blocks; ++i)
1859 for (i = 0; i < m; ++i)
1865 for (i = 1; i < blocks; ++i)
1866 beg_row[i] = beg_row[i - 1] + row_count[i - 1];
1870 printf_dh (
"\nrow_counts: ");
1871 for (i = 0; i < blocks; ++i)
1872 printf_dh (
" %i", row_count[i]);
1873 printf_dh (
"\nbeg_row: ");
1874 for (i = 0; i < blocks; ++i)
1875 printf_dh (
" %i", beg_row[i] + 1);
1881 int *tmp = (
int *) MALLOC_DH (blocks *
sizeof (
int));
1883 memcpy (tmp, beg_row, blocks *
sizeof (
int));
1884 for (i = 0; i < m; ++i)