Actual source code: matscalapack.c
1: #include <petsc/private/petscscalapack.h>
3: const char ScaLAPACKCitation[] = "@BOOK{scalapack-user-guide,\n"
4: " AUTHOR = {L. S. Blackford and J. Choi and A. Cleary and E. D'Azevedo and\n"
5: " J. Demmel and I. Dhillon and J. Dongarra and S. Hammarling and\n"
6: " G. Henry and A. Petitet and K. Stanley and D. Walker and R. C. Whaley},\n"
7: " TITLE = {Sca{LAPACK} Users' Guide},\n"
8: " PUBLISHER = {SIAM},\n"
9: " ADDRESS = {Philadelphia, PA},\n"
10: " YEAR = 1997\n"
11: "}\n";
12: static PetscBool ScaLAPACKCite = PETSC_FALSE;
14: #define DEFAULT_BLOCKSIZE 64
16: /*
17: The variable Petsc_ScaLAPACK_keyval is used to indicate an MPI attribute that
18: is attached to a communicator, in this case the attribute is a Mat_ScaLAPACK_Grid
19: */
20: static PetscMPIInt Petsc_ScaLAPACK_keyval = MPI_KEYVAL_INVALID;
22: static PetscErrorCode Petsc_ScaLAPACK_keyval_free(void)
23: {
24: PetscFunctionBegin;
25: PetscCall(PetscInfo(NULL, "Freeing Petsc_ScaLAPACK_keyval\n"));
26: PetscCallMPI(MPI_Comm_free_keyval(&Petsc_ScaLAPACK_keyval));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode MatView_ScaLAPACK(Mat A, PetscViewer viewer)
31: {
32: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
33: PetscBool isascii;
34: PetscViewerFormat format;
35: Mat Adense;
37: PetscFunctionBegin;
38: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
39: if (isascii) {
40: PetscCall(PetscViewerGetFormat(viewer, &format));
41: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
42: PetscCall(PetscViewerASCIIPrintf(viewer, "block sizes: %d,%d\n", (int)a->mb, (int)a->nb));
43: PetscCall(PetscViewerASCIIPrintf(viewer, "grid height=%d, grid width=%d\n", (int)a->grid->nprow, (int)a->grid->npcol));
44: PetscCall(PetscViewerASCIIPrintf(viewer, "coordinates of process owning first row and column: (%d,%d)\n", (int)a->rsrc, (int)a->csrc));
45: PetscCall(PetscViewerASCIIPrintf(viewer, "dimension of largest local matrix: %d x %d\n", (int)a->locr, (int)a->locc));
46: PetscFunctionReturn(PETSC_SUCCESS);
47: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
50: }
51: /* convert to dense format and call MatView() */
52: PetscCall(MatConvert(A, MATDENSE, MAT_INITIAL_MATRIX, &Adense));
53: PetscCall(MatView(Adense, viewer));
54: PetscCall(MatDestroy(&Adense));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: static PetscErrorCode MatGetInfo_ScaLAPACK(Mat A, MatInfoType flag, MatInfo *info)
59: {
60: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
61: PetscLogDouble irecv[2];
63: PetscFunctionBegin;
64: info->block_size = 1.0;
66: irecv[0] = a->lld * a->locc; /* locally allocated */
67: irecv[1] = a->locr * a->locc; /* used submatrix */
68: if (flag == MAT_LOCAL || flag == MAT_GLOBAL_MAX) {
69: info->nz_allocated = irecv[0];
70: info->nz_used = irecv[1];
71: } else if (flag == MAT_GLOBAL_MAX) {
72: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 2, MPIU_PETSCLOGDOUBLE, MPI_MAX, PetscObjectComm((PetscObject)A)));
73: info->nz_allocated = irecv[0];
74: info->nz_used = irecv[1];
75: } else if (flag == MAT_GLOBAL_SUM) {
76: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, irecv, 2, MPIU_PETSCLOGDOUBLE, MPI_SUM, PetscObjectComm((PetscObject)A)));
77: info->nz_allocated = irecv[0];
78: info->nz_used = irecv[1];
79: }
81: info->nz_unneeded = 0;
82: info->assemblies = A->num_ass;
83: info->mallocs = 0;
84: info->memory = 0; /* REVIEW ME */
85: info->fill_ratio_given = 0;
86: info->fill_ratio_needed = 0;
87: info->factor_mallocs = 0;
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: static PetscErrorCode MatSetOption_ScaLAPACK(Mat A, MatOption op, PetscBool flg)
92: {
93: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
95: PetscFunctionBegin;
96: switch (op) {
97: case MAT_NEW_NONZERO_LOCATIONS:
98: case MAT_NEW_NONZERO_LOCATION_ERR:
99: case MAT_NEW_NONZERO_ALLOCATION_ERR:
100: case MAT_SYMMETRIC:
101: case MAT_SORTED_FULL:
102: case MAT_HERMITIAN:
103: break;
104: case MAT_ROW_ORIENTED:
105: a->roworiented = flg;
106: break;
107: default:
108: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported option %s", MatOptions[op]);
109: }
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode MatSetValues_ScaLAPACK(Mat A, PetscInt nr, const PetscInt *rows, PetscInt nc, const PetscInt *cols, const PetscScalar *vals, InsertMode imode)
114: {
115: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
116: PetscInt j;
117: PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc;
118: PetscBool roworiented = a->roworiented;
120: PetscFunctionBegin;
121: PetscCheck(imode == INSERT_VALUES || imode == ADD_VALUES, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)imode);
122: for (PetscInt i = 0; i < nr; i++) {
123: if (rows[i] < 0) continue;
124: PetscCall(PetscBLASIntCast(rows[i] + 1, &gridx));
125: for (j = 0; j < nc; j++) {
126: if (cols[j] < 0) continue;
127: PetscCall(PetscBLASIntCast(cols[j] + 1, &gcidx));
128: PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc));
129: if (rsrc == a->grid->myrow && csrc == a->grid->mycol) {
130: if (roworiented) {
131: switch (imode) {
132: case INSERT_VALUES:
133: a->loc[lridx - 1 + (lcidx - 1) * a->lld] = vals[i * nc + j];
134: break;
135: default:
136: a->loc[lridx - 1 + (lcidx - 1) * a->lld] += vals[i * nc + j];
137: break;
138: }
139: } else {
140: switch (imode) {
141: case INSERT_VALUES:
142: a->loc[lridx - 1 + (lcidx - 1) * a->lld] = vals[i + j * nr];
143: break;
144: default:
145: a->loc[lridx - 1 + (lcidx - 1) * a->lld] += vals[i + j * nr];
146: break;
147: }
148: }
149: } else {
150: PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process entry even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set");
151: A->assembled = PETSC_FALSE;
152: PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], 1, cols + j, roworiented ? vals + i * nc + j : vals + i + j * nr, (PetscBool)(imode == ADD_VALUES)));
153: }
154: }
155: }
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode MatMultXXXYYY_ScaLAPACK(Mat A, PetscBool transpose, PetscBool hermitian, PetscScalar beta, const PetscScalar *x, PetscScalar *y)
160: {
161: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
162: PetscScalar *x2d, *y2d, alpha = 1.0;
163: const PetscInt *ranges;
164: PetscBLASInt xdesc[9], ydesc[9], x2desc[9], y2desc[9], mb, nb, lszx, lszy, zero = 0, one = 1, xlld, ylld;
166: PetscFunctionBegin;
167: if (transpose) {
168: /* create ScaLAPACK descriptors for vectors (1d block distribution) */
169: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
170: PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* x block size */
171: PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &xlld));
172: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(xdesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &xlld, &info));
173: PetscCall(PetscLayoutGetRanges(A->cmap, &ranges));
174: PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* y block size */
175: ylld = 1;
176: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ydesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &ylld, &info));
178: /* allocate 2d vectors */
179: lszx = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
180: lszy = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol);
181: PetscCall(PetscMalloc2(lszx, &x2d, lszy, &y2d));
182: PetscCall(PetscBLASIntCast(PetscMax(1, lszx), &xlld));
184: /* create ScaLAPACK descriptors for vectors (2d block distribution) */
185: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(x2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &xlld, &info));
186: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(y2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &ylld, &info));
188: /* redistribute x as a column of a 2d matrix */
189: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxcol));
191: /* redistribute y as a row of a 2d matrix */
192: if (beta != 0.0) PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, y, &one, &one, ydesc, y2d, &one, &one, y2desc, &a->grid->ictxrow));
194: /* call PBLAS subroutine */
195: if (hermitian) PetscCallBLAS("PBLASgemv", PBLASgemv_("C", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one));
196: else PetscCallBLAS("PBLASgemv", PBLASgemv_("T", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one));
198: /* redistribute y from a row of a 2d matrix */
199: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, y2d, &one, &one, y2desc, y, &one, &one, ydesc, &a->grid->ictxrow));
201: } else { /* non-transpose */
203: /* create ScaLAPACK descriptors for vectors (1d block distribution) */
204: PetscCall(PetscLayoutGetRanges(A->cmap, &ranges));
205: PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* x block size */
206: xlld = 1;
207: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(xdesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &xlld, &info));
208: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
209: PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* y block size */
210: PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &ylld));
211: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ydesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &ylld, &info));
213: /* allocate 2d vectors */
214: lszy = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
215: lszx = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol);
216: PetscCall(PetscMalloc2(lszx, &x2d, lszy, &y2d));
217: PetscCall(PetscBLASIntCast(PetscMax(1, lszy), &ylld));
219: /* create ScaLAPACK descriptors for vectors (2d block distribution) */
220: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(x2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &xlld, &info));
221: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(y2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &ylld, &info));
223: /* redistribute x as a row of a 2d matrix */
224: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxrow));
226: /* redistribute y as a column of a 2d matrix */
227: if (beta != 0.0) PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, y, &one, &one, ydesc, y2d, &one, &one, y2desc, &a->grid->ictxcol));
229: /* call PBLAS subroutine */
230: PetscCallBLAS("PBLASgemv", PBLASgemv_("N", &a->M, &a->N, &alpha, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &one, &beta, y2d, &one, &one, y2desc, &one));
232: /* redistribute y from a column of a 2d matrix */
233: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, y2d, &one, &one, y2desc, y, &one, &one, ydesc, &a->grid->ictxcol));
234: }
235: PetscCall(PetscFree2(x2d, y2d));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: static PetscErrorCode MatMult_ScaLAPACK(Mat A, Vec x, Vec y)
240: {
241: const PetscScalar *xarray;
242: PetscScalar *yarray;
244: PetscFunctionBegin;
245: PetscCall(VecGetArrayRead(x, &xarray));
246: PetscCall(VecGetArray(y, &yarray));
247: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_FALSE, PETSC_FALSE, 0.0, xarray, yarray));
248: PetscCall(VecRestoreArrayRead(x, &xarray));
249: PetscCall(VecRestoreArray(y, &yarray));
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: static PetscErrorCode MatMultTranspose_ScaLAPACK(Mat A, Vec x, Vec y)
254: {
255: const PetscScalar *xarray;
256: PetscScalar *yarray;
258: PetscFunctionBegin;
259: PetscCall(VecGetArrayRead(x, &xarray));
260: PetscCall(VecGetArray(y, &yarray));
261: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_FALSE, 0.0, xarray, yarray));
262: PetscCall(VecRestoreArrayRead(x, &xarray));
263: PetscCall(VecRestoreArray(y, &yarray));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode MatMultHermitianTranspose_ScaLAPACK(Mat A, Vec x, Vec y)
268: {
269: const PetscScalar *xarray;
270: PetscScalar *yarray;
272: PetscFunctionBegin;
273: PetscCall(VecGetArrayRead(x, &xarray));
274: PetscCall(VecGetArrayWrite(y, &yarray));
275: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_TRUE, 0.0, xarray, yarray));
276: PetscCall(VecRestoreArrayRead(x, &xarray));
277: PetscCall(VecRestoreArrayWrite(y, &yarray));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode MatMultAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z)
282: {
283: const PetscScalar *xarray;
284: PetscScalar *zarray;
286: PetscFunctionBegin;
287: if (y != z) PetscCall(VecCopy(y, z));
288: PetscCall(VecGetArrayRead(x, &xarray));
289: PetscCall(VecGetArray(z, &zarray));
290: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_FALSE, PETSC_FALSE, 1.0, xarray, zarray));
291: PetscCall(VecRestoreArrayRead(x, &xarray));
292: PetscCall(VecRestoreArray(z, &zarray));
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode MatMultTransposeAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z)
297: {
298: const PetscScalar *xarray;
299: PetscScalar *zarray;
301: PetscFunctionBegin;
302: if (y != z) PetscCall(VecCopy(y, z));
303: PetscCall(VecGetArrayRead(x, &xarray));
304: PetscCall(VecGetArray(z, &zarray));
305: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_FALSE, 1.0, xarray, zarray));
306: PetscCall(VecRestoreArrayRead(x, &xarray));
307: PetscCall(VecRestoreArray(z, &zarray));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: static PetscErrorCode MatMultHermitianTransposeAdd_ScaLAPACK(Mat A, Vec x, Vec y, Vec z)
312: {
313: const PetscScalar *xarray;
314: PetscScalar *zarray;
316: PetscFunctionBegin;
317: if (y != z) PetscCall(VecCopy(y, z));
318: PetscCall(VecGetArrayRead(x, &xarray));
319: PetscCall(VecGetArray(z, &zarray));
320: PetscCall(MatMultXXXYYY_ScaLAPACK(A, PETSC_TRUE, PETSC_TRUE, 1.0, xarray, zarray));
321: PetscCall(VecRestoreArrayRead(x, &xarray));
322: PetscCall(VecRestoreArray(z, &zarray));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: PetscErrorCode MatMatMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C)
327: {
328: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
329: Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data;
330: Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data;
331: PetscScalar sone = 1.0, zero = 0.0;
332: PetscBLASInt one = 1;
334: PetscFunctionBegin;
335: PetscCallBLAS("PBLASgemm", PBLASgemm_("N", "N", &a->M, &b->N, &a->N, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc));
336: C->assembled = PETSC_TRUE;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: PetscErrorCode MatMatMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C)
341: {
342: PetscFunctionBegin;
343: PetscCall(MatSetSizes(C, A->rmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
344: PetscCall(MatSetType(C, MATSCALAPACK));
345: PetscCall(MatSetUp(C));
346: C->ops->matmultnumeric = MatMatMultNumeric_ScaLAPACK;
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: static PetscErrorCode MatTransposeMatMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C)
351: {
352: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
353: Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data;
354: Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data;
355: PetscScalar sone = 1.0, zero = 0.0;
356: PetscBLASInt one = 1;
358: PetscFunctionBegin;
359: PetscCallBLAS("PBLASgemm", PBLASgemm_("T", "N", &a->N, &b->N, &a->M, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc));
360: C->assembled = PETSC_TRUE;
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: static PetscErrorCode MatTransposeMatMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C)
365: {
366: PetscFunctionBegin;
367: PetscCall(MatSetSizes(C, A->cmap->n, B->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
368: PetscCall(MatSetType(C, MATSCALAPACK));
369: PetscCall(MatSetUp(C));
370: C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_ScaLAPACK;
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static PetscErrorCode MatMatTransposeMultNumeric_ScaLAPACK(Mat A, Mat B, Mat C)
375: {
376: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
377: Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data;
378: Mat_ScaLAPACK *c = (Mat_ScaLAPACK *)C->data;
379: PetscScalar sone = 1.0, zero = 0.0;
380: PetscBLASInt one = 1;
382: PetscFunctionBegin;
383: PetscCallBLAS("PBLASgemm", PBLASgemm_("N", "T", &a->M, &b->M, &a->N, &sone, a->loc, &one, &one, a->desc, b->loc, &one, &one, b->desc, &zero, c->loc, &one, &one, c->desc));
384: C->assembled = PETSC_TRUE;
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode MatMatTransposeMultSymbolic_ScaLAPACK(Mat A, Mat B, PetscReal fill, Mat C)
389: {
390: PetscFunctionBegin;
391: PetscCall(MatSetSizes(C, A->rmap->n, B->rmap->n, PETSC_DECIDE, PETSC_DECIDE));
392: PetscCall(MatSetType(C, MATSCALAPACK));
393: PetscCall(MatSetUp(C));
394: C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_ScaLAPACK;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_ScaLAPACK(Mat C)
399: {
400: Mat_Product *product = C->product;
402: PetscFunctionBegin;
403: switch (product->type) {
404: case MATPRODUCT_AB:
405: C->ops->matmultsymbolic = MatMatMultSymbolic_ScaLAPACK;
406: C->ops->productsymbolic = MatProductSymbolic_AB;
407: break;
408: case MATPRODUCT_AtB:
409: C->ops->transposematmultsymbolic = MatTransposeMatMultSymbolic_ScaLAPACK;
410: C->ops->productsymbolic = MatProductSymbolic_AtB;
411: break;
412: case MATPRODUCT_ABt:
413: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_ScaLAPACK;
414: C->ops->productsymbolic = MatProductSymbolic_ABt;
415: break;
416: default:
417: SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatProduct type %s is not supported for ScaLAPACK and ScaLAPACK matrices", MatProductTypes[product->type]);
418: }
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode MatGetDiagonal_ScaLAPACK(Mat A, Vec D)
423: {
424: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
425: PetscScalar *darray, *d2d, v;
426: const PetscInt *ranges;
427: PetscBLASInt j, ddesc[9], d2desc[9], mb, nb, lszd, zero = 0, one = 1, dlld;
429: PetscFunctionBegin;
430: PetscCall(VecGetArray(D, &darray));
432: if (A->rmap->N <= A->cmap->N) { /* row version */
434: /* create ScaLAPACK descriptor for vector (1d block distribution) */
435: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
436: PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* D block size */
437: PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &dlld));
438: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ddesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &dlld, &info));
440: /* allocate 2d vector */
441: lszd = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
442: PetscCall(PetscCalloc1(lszd, &d2d));
443: PetscCall(PetscBLASIntCast(PetscMax(1, lszd), &dlld));
445: /* create ScaLAPACK descriptor for vector (2d block distribution) */
446: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(d2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &dlld, &info));
448: /* collect diagonal */
449: for (j = 1; j <= a->M; j++) {
450: PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("R", " ", &v, a->loc, &j, &j, a->desc));
451: PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(d2d, &j, &one, d2desc, &v));
452: }
454: /* redistribute d from a column of a 2d matrix */
455: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, d2d, &one, &one, d2desc, darray, &one, &one, ddesc, &a->grid->ictxcol));
456: PetscCall(PetscFree(d2d));
458: } else { /* column version */
460: /* create ScaLAPACK descriptor for vector (1d block distribution) */
461: PetscCall(PetscLayoutGetRanges(A->cmap, &ranges));
462: PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* D block size */
463: dlld = 1;
464: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ddesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &dlld, &info));
466: /* allocate 2d vector */
467: lszd = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol);
468: PetscCall(PetscCalloc1(lszd, &d2d));
470: /* create ScaLAPACK descriptor for vector (2d block distribution) */
471: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(d2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &dlld, &info));
473: /* collect diagonal */
474: for (j = 1; j <= a->N; j++) {
475: PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("C", " ", &v, a->loc, &j, &j, a->desc));
476: PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(d2d, &one, &j, d2desc, &v));
477: }
479: /* redistribute d from a row of a 2d matrix */
480: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, d2d, &one, &one, d2desc, darray, &one, &one, ddesc, &a->grid->ictxrow));
481: PetscCall(PetscFree(d2d));
482: }
484: PetscCall(VecRestoreArray(D, &darray));
485: PetscCall(VecAssemblyBegin(D));
486: PetscCall(VecAssemblyEnd(D));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode MatDiagonalScale_ScaLAPACK(Mat A, Vec L, Vec R)
491: {
492: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
493: const PetscScalar *d;
494: const PetscInt *ranges;
495: PetscScalar *d2d;
496: PetscBLASInt i, j, ddesc[9], d2desc[9], mb, nb, lszd, zero = 0, one = 1, dlld;
498: PetscFunctionBegin;
499: if (R) {
500: PetscCall(VecGetArrayRead(R, &d));
501: /* create ScaLAPACK descriptor for vector (1d block distribution) */
502: PetscCall(PetscLayoutGetRanges(A->cmap, &ranges));
503: PetscCall(PetscBLASIntCast(ranges[1], &nb)); /* D block size */
504: dlld = 1;
505: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ddesc, &one, &a->N, &one, &nb, &zero, &zero, &a->grid->ictxrow, &dlld, &info));
507: /* allocate 2d vector */
508: lszd = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol);
509: PetscCall(PetscCalloc1(lszd, &d2d));
511: /* create ScaLAPACK descriptor for vector (2d block distribution) */
512: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(d2desc, &one, &a->N, &one, &a->nb, &zero, &zero, &a->grid->ictxt, &dlld, &info));
514: /* redistribute d to a row of a 2d matrix */
515: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&one, &a->N, d, &one, &one, ddesc, d2d, &one, &one, d2desc, &a->grid->ictxrow));
517: /* broadcast along process columns */
518: if (!a->grid->myrow) Cdgebs2d(a->grid->ictxt, "C", " ", 1, lszd, d2d, dlld);
519: else Cdgebr2d(a->grid->ictxt, "C", " ", 1, lszd, d2d, dlld, 0, a->grid->mycol);
521: /* local scaling */
522: for (j = 0; j < a->locc; j++)
523: for (i = 0; i < a->locr; i++) a->loc[i + j * a->lld] *= d2d[j];
525: PetscCall(PetscFree(d2d));
526: PetscCall(VecRestoreArrayRead(R, &d));
527: }
528: if (L) {
529: PetscCall(VecGetArrayRead(L, &d));
530: /* create ScaLAPACK descriptor for vector (1d block distribution) */
531: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
532: PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* D block size */
533: PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &dlld));
534: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(ddesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &dlld, &info));
536: /* allocate 2d vector */
537: lszd = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
538: PetscCall(PetscCalloc1(lszd, &d2d));
539: PetscCall(PetscBLASIntCast(PetscMax(1, lszd), &dlld));
541: /* create ScaLAPACK descriptor for vector (2d block distribution) */
542: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(d2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &dlld, &info));
544: /* redistribute d to a column of a 2d matrix */
545: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, d, &one, &one, ddesc, d2d, &one, &one, d2desc, &a->grid->ictxcol));
547: /* broadcast along process rows */
548: if (!a->grid->mycol) Cdgebs2d(a->grid->ictxt, "R", " ", lszd, 1, d2d, dlld);
549: else Cdgebr2d(a->grid->ictxt, "R", " ", lszd, 1, d2d, dlld, a->grid->myrow, 0);
551: /* local scaling */
552: for (i = 0; i < a->locr; i++)
553: for (j = 0; j < a->locc; j++) a->loc[i + j * a->lld] *= d2d[i];
555: PetscCall(PetscFree(d2d));
556: PetscCall(VecRestoreArrayRead(L, &d));
557: }
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: static PetscErrorCode MatScale_ScaLAPACK(Mat X, PetscScalar a)
562: {
563: Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data;
564: PetscBLASInt n, one = 1;
566: PetscFunctionBegin;
567: n = x->lld * x->locc;
568: PetscCallBLAS("BLASscal", BLASscal_(&n, &a, x->loc, &one));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode MatShift_ScaLAPACK(Mat X, PetscScalar alpha)
573: {
574: Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data;
575: PetscBLASInt i, n;
576: PetscScalar v;
578: PetscFunctionBegin;
579: n = PetscMin(x->M, x->N);
580: for (i = 1; i <= n; i++) {
581: PetscCallBLAS("SCALAPACKelget", SCALAPACKelget_("-", " ", &v, x->loc, &i, &i, x->desc));
582: v += alpha;
583: PetscCallBLAS("SCALAPACKelset", SCALAPACKelset_(x->loc, &i, &i, x->desc, &v));
584: }
585: PetscFunctionReturn(PETSC_SUCCESS);
586: }
588: static PetscErrorCode MatAXPY_ScaLAPACK(Mat Y, PetscScalar alpha, Mat X, MatStructure str)
589: {
590: Mat_ScaLAPACK *x = (Mat_ScaLAPACK *)X->data;
591: Mat_ScaLAPACK *y = (Mat_ScaLAPACK *)Y->data;
592: PetscBLASInt one = 1;
593: PetscScalar beta = 1.0;
595: PetscFunctionBegin;
596: MatScaLAPACKCheckDistribution(Y, 1, X, 3);
597: PetscCallBLAS("SCALAPACKmatadd", SCALAPACKmatadd_(&x->M, &x->N, &alpha, x->loc, &one, &one, x->desc, &beta, y->loc, &one, &one, y->desc));
598: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode MatCopy_ScaLAPACK(Mat A, Mat B, MatStructure str)
603: {
604: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
605: Mat_ScaLAPACK *b = (Mat_ScaLAPACK *)B->data;
607: PetscFunctionBegin;
608: PetscCall(PetscArraycpy(b->loc, a->loc, a->lld * a->locc));
609: PetscCall(PetscObjectStateIncrease((PetscObject)B));
610: PetscFunctionReturn(PETSC_SUCCESS);
611: }
613: static PetscErrorCode MatDuplicate_ScaLAPACK(Mat A, MatDuplicateOption op, Mat *B)
614: {
615: Mat Bs;
616: MPI_Comm comm;
617: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b;
619: PetscFunctionBegin;
620: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
621: PetscCall(MatCreate(comm, &Bs));
622: PetscCall(MatSetSizes(Bs, A->rmap->n, A->cmap->n, PETSC_DECIDE, PETSC_DECIDE));
623: PetscCall(MatSetType(Bs, MATSCALAPACK));
624: b = (Mat_ScaLAPACK *)Bs->data;
625: b->M = a->M;
626: b->N = a->N;
627: b->mb = a->mb;
628: b->nb = a->nb;
629: b->rsrc = a->rsrc;
630: b->csrc = a->csrc;
631: PetscCall(MatSetUp(Bs));
632: *B = Bs;
633: if (op == MAT_COPY_VALUES) PetscCall(PetscArraycpy(b->loc, a->loc, a->lld * a->locc));
634: Bs->assembled = PETSC_TRUE;
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: static PetscErrorCode MatTranspose_ScaLAPACK(Mat A, MatReuse reuse, Mat *B)
639: {
640: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b;
641: Mat Bs = *B;
642: PetscBLASInt one = 1;
643: PetscScalar sone = 1.0, zero = 0.0;
644: #if PetscDefined(USE_COMPLEX)
645: PetscInt i;
646: #endif
648: PetscFunctionBegin;
649: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
650: PetscCheck(reuse == MAT_INITIAL_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_INITIAL_MATRIX supported");
651: PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->nb, a->mb, a->N, a->M, a->csrc, a->rsrc, &Bs));
652: *B = Bs;
653: b = (Mat_ScaLAPACK *)Bs->data;
654: PetscCallBLAS("PBLAStran", PBLAStran_(&a->N, &a->M, &sone, a->loc, &one, &one, a->desc, &zero, b->loc, &one, &one, b->desc));
655: #if PetscDefined(USE_COMPLEX)
656: /* undo conjugation */
657: for (i = 0; i < b->locr; i++)
658: for (PetscInt j = 0; j < b->locc; j++) b->loc[i + j * b->lld] = PetscConj(b->loc[i + j * b->lld]);
659: #endif
660: Bs->assembled = PETSC_TRUE;
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: static PetscErrorCode MatConjugate_ScaLAPACK(Mat A)
665: {
666: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
667: PetscInt i;
669: PetscFunctionBegin;
670: for (i = 0; i < a->locr; i++)
671: for (PetscInt j = 0; j < a->locc; j++) a->loc[i + j * a->lld] = PetscConj(a->loc[i + j * a->lld]);
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: static PetscErrorCode MatHermitianTranspose_ScaLAPACK(Mat A, MatReuse reuse, Mat *B)
676: {
677: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *b;
678: Mat Bs = *B;
679: PetscBLASInt one = 1;
680: PetscScalar sone = 1.0, zero = 0.0;
682: PetscFunctionBegin;
683: PetscCheck(reuse == MAT_INITIAL_MATRIX, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only MAT_INITIAL_MATRIX supported");
684: PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->nb, a->mb, a->N, a->M, a->csrc, a->rsrc, &Bs));
685: *B = Bs;
686: b = (Mat_ScaLAPACK *)Bs->data;
687: PetscCallBLAS("PBLAStran", PBLAStran_(&a->N, &a->M, &sone, a->loc, &one, &one, a->desc, &zero, b->loc, &one, &one, b->desc));
688: Bs->assembled = PETSC_TRUE;
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: static PetscErrorCode MatSolve_ScaLAPACK(Mat A, Vec B, Vec X)
693: {
694: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
695: PetscScalar *x, *x2d;
696: const PetscInt *ranges;
697: PetscBLASInt xdesc[9], x2desc[9], mb, lszx, zero = 0, one = 1, xlld, nrhs = 1;
699: PetscFunctionBegin;
700: PetscCall(VecCopy(B, X));
701: PetscCall(VecGetArray(X, &x));
703: /* create ScaLAPACK descriptor for a vector (1d block distribution) */
704: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
705: PetscCall(PetscBLASIntCast(ranges[1], &mb)); /* x block size */
706: PetscCall(PetscBLASIntCast(PetscMax(1, A->rmap->n), &xlld));
707: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(xdesc, &a->M, &one, &mb, &one, &zero, &zero, &a->grid->ictxcol, &xlld, &info));
709: /* allocate 2d vector */
710: lszx = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
711: PetscCall(PetscMalloc1(lszx, &x2d));
712: PetscCall(PetscBLASIntCast(PetscMax(1, lszx), &xlld));
714: /* create ScaLAPACK descriptor for a vector (2d block distribution) */
715: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(x2desc, &a->M, &one, &a->mb, &one, &zero, &zero, &a->grid->ictxt, &xlld, &info));
717: /* redistribute x as a column of a 2d matrix */
718: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x, &one, &one, xdesc, x2d, &one, &one, x2desc, &a->grid->ictxcol));
720: /* call ScaLAPACK subroutine */
721: switch (A->factortype) {
722: case MAT_FACTOR_LU:
723: PetscCallScaLAPACKInfo("getrs", SCALAPACKgetrs_("N", &a->M, &nrhs, a->loc, &one, &one, a->desc, a->pivots, x2d, &one, &one, x2desc, &info));
724: break;
725: case MAT_FACTOR_CHOLESKY:
726: PetscCallScaLAPACKInfo("potrs", SCALAPACKpotrs_("L", &a->M, &nrhs, a->loc, &one, &one, a->desc, x2d, &one, &one, x2desc, &info));
727: break;
728: default:
729: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
730: }
732: /* redistribute x from a column of a 2d matrix */
733: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &one, x2d, &one, &one, x2desc, x, &one, &one, xdesc, &a->grid->ictxcol));
735: PetscCall(PetscFree(x2d));
736: PetscCall(VecRestoreArray(X, &x));
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: static PetscErrorCode MatSolveAdd_ScaLAPACK(Mat A, Vec B, Vec Y, Vec X)
741: {
742: PetscFunctionBegin;
743: PetscCall(MatSolve_ScaLAPACK(A, B, X));
744: PetscCall(VecAXPY(X, 1, Y));
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: static PetscErrorCode MatMatSolve_ScaLAPACK(Mat A, Mat B, Mat X)
749: {
750: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data, *x;
751: PetscBool flg1, flg2;
752: PetscBLASInt one = 1;
753: Mat C;
754: MatType type;
756: PetscFunctionBegin;
757: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSCALAPACK, &flg1));
758: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATSCALAPACK, &flg2));
759: if (flg1 && flg2) MatScaLAPACKCheckDistribution(B, 2, X, 3);
760: if (flg2) {
761: if (flg1) PetscCall(MatCopy(B, X, SAME_NONZERO_PATTERN));
762: else PetscCall(MatConvert(B, MATSCALAPACK, MAT_REUSE_MATRIX, &X));
763: C = X;
764: } else {
765: PetscCall(MatConvert(B, MATSCALAPACK, MAT_INITIAL_MATRIX, &C));
766: }
767: x = (Mat_ScaLAPACK *)C->data;
769: switch (A->factortype) {
770: case MAT_FACTOR_LU:
771: PetscCallScaLAPACKInfo("getrs", SCALAPACKgetrs_("N", &a->M, &x->N, a->loc, &one, &one, a->desc, a->pivots, x->loc, &one, &one, x->desc, &info));
772: break;
773: case MAT_FACTOR_CHOLESKY:
774: PetscCallScaLAPACKInfo("potrs", SCALAPACKpotrs_("L", &a->M, &x->N, a->loc, &one, &one, a->desc, x->loc, &one, &one, x->desc, &info));
775: break;
776: default:
777: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unfactored Matrix or Unsupported MatFactorType");
778: }
779: if (!flg2) {
780: PetscCall(MatGetType(X, &type));
781: PetscCall(MatConvert(C, type, MAT_REUSE_MATRIX, &X));
782: PetscCall(MatDestroy(&C));
783: }
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }
787: static PetscErrorCode MatLUFactor_ScaLAPACK(Mat A, IS row, IS col, const MatFactorInfo *factorinfo)
788: {
789: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
790: PetscBLASInt one = 1;
792: PetscFunctionBegin;
793: if (!a->pivots) PetscCall(PetscMalloc1(a->locr + a->mb, &a->pivots));
794: PetscCallScaLAPACKInfo("getrf", SCALAPACKgetrf_(&a->M, &a->N, a->loc, &one, &one, a->desc, a->pivots, &info));
795: A->factortype = MAT_FACTOR_LU;
796: A->assembled = PETSC_TRUE;
798: PetscCall(PetscFree(A->solvertype));
799: PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &A->solvertype));
800: PetscFunctionReturn(PETSC_SUCCESS);
801: }
803: static PetscErrorCode MatLUFactorNumeric_ScaLAPACK(Mat F, Mat A, const MatFactorInfo *info)
804: {
805: PetscFunctionBegin;
806: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
807: PetscCall(MatLUFactor_ScaLAPACK(F, 0, 0, info));
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
811: static PetscErrorCode MatLUFactorSymbolic_ScaLAPACK(Mat F, Mat A, IS r, IS c, const MatFactorInfo *info)
812: {
813: PetscFunctionBegin;
814: /* F is created and allocated by MatGetFactor_scalapack_petsc(), skip this routine. */
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: static PetscErrorCode MatCholeskyFactor_ScaLAPACK(Mat A, IS perm, const MatFactorInfo *factorinfo)
819: {
820: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
821: PetscBLASInt one = 1;
823: PetscFunctionBegin;
824: PetscCallScaLAPACKInfo("potrf", SCALAPACKpotrf_("L", &a->M, a->loc, &one, &one, a->desc, &info));
825: A->factortype = MAT_FACTOR_CHOLESKY;
826: A->assembled = PETSC_TRUE;
828: PetscCall(PetscFree(A->solvertype));
829: PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &A->solvertype));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: static PetscErrorCode MatCholeskyFactorNumeric_ScaLAPACK(Mat F, Mat A, const MatFactorInfo *info)
834: {
835: PetscFunctionBegin;
836: PetscCall(MatCopy(A, F, SAME_NONZERO_PATTERN));
837: PetscCall(MatCholeskyFactor_ScaLAPACK(F, 0, info));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode MatCholeskyFactorSymbolic_ScaLAPACK(Mat F, Mat A, IS perm, const MatFactorInfo *info)
842: {
843: PetscFunctionBegin;
844: /* F is created and allocated by MatGetFactor_scalapack_petsc(), skip this routine. */
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static PetscErrorCode MatFactorGetSolverType_scalapack_scalapack(Mat A, MatSolverType *type)
849: {
850: PetscFunctionBegin;
851: *type = MATSOLVERSCALAPACK;
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode MatGetFactor_scalapack_scalapack(Mat A, MatFactorType ftype, Mat *F)
856: {
857: Mat B;
858: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
860: PetscFunctionBegin;
861: /* Create the factorization matrix */
862: PetscCall(MatCreateScaLAPACK(PetscObjectComm((PetscObject)A), a->mb, a->nb, a->M, a->N, a->rsrc, a->csrc, &B));
863: B->trivialsymbolic = PETSC_TRUE;
864: B->factortype = ftype;
865: PetscCall(PetscFree(B->solvertype));
866: PetscCall(PetscStrallocpy(MATSOLVERSCALAPACK, &B->solvertype));
868: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatFactorGetSolverType_C", MatFactorGetSolverType_scalapack_scalapack));
869: *F = B;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: PETSC_INTERN PetscErrorCode MatSolverTypeRegister_ScaLAPACK(void)
874: {
875: PetscFunctionBegin;
876: PetscCall(MatSolverTypeRegister(MATSOLVERSCALAPACK, MATSCALAPACK, MAT_FACTOR_LU, MatGetFactor_scalapack_scalapack));
877: PetscCall(MatSolverTypeRegister(MATSOLVERSCALAPACK, MATSCALAPACK, MAT_FACTOR_CHOLESKY, MatGetFactor_scalapack_scalapack));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode MatNorm_ScaLAPACK(Mat A, NormType type, PetscReal *nrm)
882: {
883: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
884: PetscBLASInt one = 1, lwork = 0;
885: const char *ntype;
886: PetscScalar *work = NULL, dummy;
888: PetscFunctionBegin;
889: switch (type) {
890: case NORM_1:
891: ntype = "1";
892: lwork = PetscMax(a->locr, a->locc);
893: break;
894: case NORM_FROBENIUS:
895: ntype = "F";
896: work = &dummy;
897: break;
898: case NORM_INFINITY:
899: ntype = "I";
900: lwork = PetscMax(a->locr, a->locc);
901: break;
902: default:
903: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Unsupported norm type");
904: }
905: if (lwork) PetscCall(PetscMalloc1(lwork, &work));
906: *nrm = SCALAPACKlange_(ntype, &a->M, &a->N, a->loc, &one, &one, a->desc, work);
907: if (lwork) PetscCall(PetscFree(work));
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: static PetscErrorCode MatZeroEntries_ScaLAPACK(Mat A)
912: {
913: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
915: PetscFunctionBegin;
916: PetscCall(PetscArrayzero(a->loc, a->lld * a->locc));
917: PetscFunctionReturn(PETSC_SUCCESS);
918: }
920: static PetscErrorCode MatGetOwnershipIS_ScaLAPACK(Mat A, IS *rows, IS *cols)
921: {
922: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
923: PetscInt i, n, nb, isrc, nproc, iproc, *idx;
925: PetscFunctionBegin;
926: if (rows) {
927: n = a->locr;
928: nb = a->mb;
929: isrc = a->rsrc;
930: nproc = a->grid->nprow;
931: iproc = a->grid->myrow;
932: PetscCall(PetscMalloc1(n, &idx));
933: for (i = 0; i < n; i++) idx[i] = nproc * nb * (i / nb) + i % nb + ((nproc + iproc - isrc) % nproc) * nb;
934: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_OWN_POINTER, rows));
935: }
936: if (cols) {
937: n = a->locc;
938: nb = a->nb;
939: isrc = a->csrc;
940: nproc = a->grid->npcol;
941: iproc = a->grid->mycol;
942: PetscCall(PetscMalloc1(n, &idx));
943: for (i = 0; i < n; i++) idx[i] = nproc * nb * (i / nb) + i % nb + ((nproc + iproc - isrc) % nproc) * nb;
944: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_OWN_POINTER, cols));
945: }
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: static PetscErrorCode MatConvert_ScaLAPACK_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *B)
950: {
951: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
952: Mat Bmpi;
953: MPI_Comm comm;
954: PetscInt i, M = A->rmap->N, N = A->cmap->N, m, n, rstart, rend, nz, ldb;
955: const PetscInt *ranges, *branges, *cwork;
956: const PetscScalar *vwork;
957: PetscBLASInt bdesc[9], bmb, zero = 0, one = 1, lld;
958: PetscScalar *barray;
959: PetscBool differ = PETSC_FALSE;
960: PetscMPIInt size;
962: PetscFunctionBegin;
963: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
964: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
966: if (reuse == MAT_REUSE_MATRIX) { /* check if local sizes differ in A and B */
967: PetscCallMPI(MPI_Comm_size(comm, &size));
968: PetscCall(PetscLayoutGetRanges((*B)->rmap, &branges));
969: for (i = 0; i < size; i++)
970: if (ranges[i + 1] != branges[i + 1]) {
971: differ = PETSC_TRUE;
972: break;
973: }
974: }
976: if (reuse == MAT_REUSE_MATRIX && differ) { /* special case, use auxiliary dense matrix */
977: PetscCall(MatCreate(comm, &Bmpi));
978: m = PETSC_DECIDE;
979: PetscCall(PetscSplitOwnershipEqual(comm, &m, &M));
980: n = PETSC_DECIDE;
981: PetscCall(PetscSplitOwnershipEqual(comm, &n, &N));
982: PetscCall(MatSetSizes(Bmpi, m, n, M, N));
983: PetscCall(MatSetType(Bmpi, MATDENSE));
984: PetscCall(MatSetUp(Bmpi));
986: /* create ScaLAPACK descriptor for B (1d block distribution) */
987: PetscCall(PetscBLASIntCast(ranges[1], &bmb)); /* row block size */
988: PetscCall(MatDenseGetLDA(Bmpi, &ldb));
989: PetscCall(PetscBLASIntCast(PetscMax(ldb, 1), &lld)); /* local leading dimension */
990: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(bdesc, &a->M, &a->N, &bmb, &a->N, &zero, &zero, &a->grid->ictxcol, &lld, &info));
992: /* redistribute matrix */
993: PetscCall(MatDenseGetArray(Bmpi, &barray));
994: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &a->N, a->loc, &one, &one, a->desc, barray, &one, &one, bdesc, &a->grid->ictxcol));
995: PetscCall(MatDenseRestoreArray(Bmpi, &barray));
997: /* transfer rows of auxiliary matrix to the final matrix B */
998: PetscCall(MatGetOwnershipRange(Bmpi, &rstart, &rend));
999: for (i = rstart; i < rend; i++) {
1000: PetscCall(MatGetRow(Bmpi, i, &nz, &cwork, &vwork));
1001: PetscCall(MatSetValues(*B, 1, &i, nz, cwork, vwork, INSERT_VALUES));
1002: PetscCall(MatRestoreRow(Bmpi, i, &nz, &cwork, &vwork));
1003: }
1004: PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1005: PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1006: PetscCall(MatDestroy(&Bmpi));
1008: } else { /* normal cases */
1010: if (reuse == MAT_REUSE_MATRIX) Bmpi = *B;
1011: else {
1012: PetscCall(MatCreate(comm, &Bmpi));
1013: m = PETSC_DECIDE;
1014: PetscCall(PetscSplitOwnershipEqual(comm, &m, &M));
1015: n = PETSC_DECIDE;
1016: PetscCall(PetscSplitOwnershipEqual(comm, &n, &N));
1017: PetscCall(MatSetSizes(Bmpi, m, n, M, N));
1018: PetscCall(MatSetType(Bmpi, MATDENSE));
1019: PetscCall(MatSetUp(Bmpi));
1020: }
1022: /* create ScaLAPACK descriptor for B (1d block distribution) */
1023: PetscCall(PetscBLASIntCast(ranges[1], &bmb)); /* row block size */
1024: PetscCall(MatDenseGetLDA(Bmpi, &ldb));
1025: PetscCall(PetscBLASIntCast(PetscMax(ldb, 1), &lld)); /* local leading dimension */
1026: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(bdesc, &a->M, &a->N, &bmb, &a->N, &zero, &zero, &a->grid->ictxcol, &lld, &info));
1028: /* redistribute matrix */
1029: PetscCall(MatDenseGetArray(Bmpi, &barray));
1030: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&a->M, &a->N, a->loc, &one, &one, a->desc, barray, &one, &one, bdesc, &a->grid->ictxcol));
1031: PetscCall(MatDenseRestoreArray(Bmpi, &barray));
1033: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &Bmpi));
1034: else *B = Bmpi;
1035: }
1036: PetscFunctionReturn(PETSC_SUCCESS);
1037: }
1039: static inline PetscErrorCode MatScaLAPACKCheckLayout(PetscLayout map, PetscBool *correct)
1040: {
1041: const PetscInt *ranges;
1042: PetscMPIInt size;
1043: PetscInt i, n;
1045: PetscFunctionBegin;
1046: *correct = PETSC_TRUE;
1047: PetscCallMPI(MPI_Comm_size(map->comm, &size));
1048: if (size > 1) {
1049: PetscCall(PetscLayoutGetRanges(map, &ranges));
1050: n = ranges[1] - ranges[0];
1051: for (i = 1; i < size; i++)
1052: if (ranges[i + 1] - ranges[i] != n) break;
1053: *correct = (PetscBool)(i == size || (i == size - 1 && ranges[i + 1] - ranges[i] <= n));
1054: }
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: PETSC_INTERN PetscErrorCode MatConvert_Dense_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *B)
1059: {
1060: Mat_ScaLAPACK *b;
1061: Mat Bmpi;
1062: MPI_Comm comm;
1063: PetscInt M = A->rmap->N, N = A->cmap->N, m, n;
1064: const PetscInt *ranges, *rows, *cols;
1065: PetscBLASInt adesc[9], amb, zero = 0, one = 1, lld;
1066: const PetscScalar *aarray;
1067: IS ir, ic;
1068: PetscInt lda;
1069: PetscBool flg;
1071: PetscFunctionBegin;
1072: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
1074: if (reuse == MAT_REUSE_MATRIX) Bmpi = *B;
1075: else {
1076: PetscCall(MatCreate(comm, &Bmpi));
1077: m = PETSC_DECIDE;
1078: PetscCall(PetscSplitOwnershipEqual(comm, &m, &M));
1079: n = PETSC_DECIDE;
1080: PetscCall(PetscSplitOwnershipEqual(comm, &n, &N));
1081: PetscCall(MatSetSizes(Bmpi, m, n, M, N));
1082: PetscCall(MatSetType(Bmpi, MATSCALAPACK));
1083: PetscCall(MatSetUp(Bmpi));
1084: }
1085: b = (Mat_ScaLAPACK *)Bmpi->data;
1087: PetscCall(MatDenseGetLDA(A, &lda));
1088: PetscCall(MatDenseGetArrayRead(A, &aarray));
1089: PetscCall(MatScaLAPACKCheckLayout(A->rmap, &flg));
1090: if (flg) PetscCall(MatScaLAPACKCheckLayout(A->cmap, &flg));
1091: if (flg) { /* if the input Mat has a ScaLAPACK-compatible layout, use ScaLAPACK for the redistribution */
1092: /* create ScaLAPACK descriptor for A (1d block distribution) */
1093: PetscCall(PetscLayoutGetRanges(A->rmap, &ranges));
1094: PetscCall(PetscBLASIntCast(ranges[1], &amb)); /* row block size */
1095: PetscCall(PetscBLASIntCast(PetscMax(lda, 1), &lld)); /* local leading dimension */
1096: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(adesc, &b->M, &b->N, &amb, &b->N, &zero, &zero, &b->grid->ictxcol, &lld, &info));
1098: /* redistribute matrix */
1099: PetscCallBLAS("SCALAPACKgemr2d", SCALAPACKgemr2d_(&b->M, &b->N, aarray, &one, &one, adesc, b->loc, &one, &one, b->desc, &b->grid->ictxcol));
1100: Bmpi->nooffprocentries = PETSC_TRUE;
1101: } else { /* if the input Mat has a ScaLAPACK-incompatible layout, redistribute via MatSetValues() */
1102: PetscCheck(lda == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Leading dimension (%" PetscInt_FMT ") different than local number of rows (%" PetscInt_FMT ")", lda, A->rmap->n);
1103: b->roworiented = PETSC_FALSE;
1104: PetscCall(MatGetOwnershipIS(A, &ir, &ic));
1105: PetscCall(ISGetIndices(ir, &rows));
1106: PetscCall(ISGetIndices(ic, &cols));
1107: PetscCall(MatSetValues(Bmpi, A->rmap->n, rows, A->cmap->N, cols, aarray, INSERT_VALUES));
1108: PetscCall(ISRestoreIndices(ir, &rows));
1109: PetscCall(ISRestoreIndices(ic, &cols));
1110: PetscCall(ISDestroy(&ic));
1111: PetscCall(ISDestroy(&ir));
1112: }
1113: PetscCall(MatDenseRestoreArrayRead(A, &aarray));
1114: PetscCall(MatAssemblyBegin(Bmpi, MAT_FINAL_ASSEMBLY));
1115: PetscCall(MatAssemblyEnd(Bmpi, MAT_FINAL_ASSEMBLY));
1116: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &Bmpi));
1117: else *B = Bmpi;
1118: PetscFunctionReturn(PETSC_SUCCESS);
1119: }
1121: PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1122: {
1123: Mat mat_scal;
1124: PetscInt M = A->rmap->N, N = A->cmap->N, rstart = A->rmap->rstart, rend = A->rmap->rend, m, n, row, ncols;
1125: const PetscInt *cols;
1126: const PetscScalar *vals;
1128: PetscFunctionBegin;
1129: if (reuse == MAT_REUSE_MATRIX) {
1130: mat_scal = *newmat;
1131: PetscCall(MatZeroEntries(mat_scal));
1132: } else {
1133: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_scal));
1134: m = PETSC_DECIDE;
1135: PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &m, &M));
1136: n = PETSC_DECIDE;
1137: PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &n, &N));
1138: PetscCall(MatSetSizes(mat_scal, m, n, M, N));
1139: PetscCall(MatSetType(mat_scal, MATSCALAPACK));
1140: PetscCall(MatSetUp(mat_scal));
1141: }
1142: for (row = rstart; row < rend; row++) {
1143: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1144: PetscCall(MatSetValues(mat_scal, 1, &row, ncols, cols, vals, INSERT_VALUES));
1145: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1146: }
1147: PetscCall(MatAssemblyBegin(mat_scal, MAT_FINAL_ASSEMBLY));
1148: PetscCall(MatAssemblyEnd(mat_scal, MAT_FINAL_ASSEMBLY));
1150: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &mat_scal));
1151: else *newmat = mat_scal;
1152: PetscFunctionReturn(PETSC_SUCCESS);
1153: }
1155: PETSC_INTERN PetscErrorCode MatConvert_SBAIJ_ScaLAPACK(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1156: {
1157: Mat mat_scal;
1158: PetscInt M = A->rmap->N, N = A->cmap->N, m, n, row, ncols, j, rstart = A->rmap->rstart, rend = A->rmap->rend;
1159: const PetscInt *cols;
1160: const PetscScalar *vals;
1161: PetscScalar v;
1163: PetscFunctionBegin;
1164: if (reuse == MAT_REUSE_MATRIX) {
1165: mat_scal = *newmat;
1166: PetscCall(MatZeroEntries(mat_scal));
1167: } else {
1168: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &mat_scal));
1169: m = PETSC_DECIDE;
1170: PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &m, &M));
1171: n = PETSC_DECIDE;
1172: PetscCall(PetscSplitOwnershipEqual(PetscObjectComm((PetscObject)A), &n, &N));
1173: PetscCall(MatSetSizes(mat_scal, m, n, M, N));
1174: PetscCall(MatSetType(mat_scal, MATSCALAPACK));
1175: PetscCall(MatSetUp(mat_scal));
1176: }
1177: PetscCall(MatGetRowUpperTriangular(A));
1178: for (row = rstart; row < rend; row++) {
1179: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
1180: PetscCall(MatSetValues(mat_scal, 1, &row, ncols, cols, vals, ADD_VALUES));
1181: for (j = 0; j < ncols; j++) { /* lower triangular part */
1182: if (cols[j] == row) continue;
1183: v = A->hermitian == PETSC_BOOL3_TRUE ? PetscConj(vals[j]) : vals[j];
1184: PetscCall(MatSetValues(mat_scal, 1, &cols[j], 1, &row, &v, ADD_VALUES));
1185: }
1186: PetscCall(MatRestoreRow(A, row, &ncols, &cols, &vals));
1187: }
1188: PetscCall(MatRestoreRowUpperTriangular(A));
1189: PetscCall(MatAssemblyBegin(mat_scal, MAT_FINAL_ASSEMBLY));
1190: PetscCall(MatAssemblyEnd(mat_scal, MAT_FINAL_ASSEMBLY));
1192: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &mat_scal));
1193: else *newmat = mat_scal;
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: static PetscErrorCode MatScaLAPACKSetPreallocation(Mat A)
1198: {
1199: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1200: PetscInt sz = 0;
1202: PetscFunctionBegin;
1203: PetscCall(PetscLayoutSetUp(A->rmap));
1204: PetscCall(PetscLayoutSetUp(A->cmap));
1205: if (!a->lld) a->lld = a->locr;
1207: PetscCall(PetscFree(a->loc));
1208: PetscCall(PetscIntMultError(a->lld, a->locc, &sz));
1209: PetscCall(PetscCalloc1(sz, &a->loc));
1211: A->preallocated = PETSC_TRUE;
1212: PetscFunctionReturn(PETSC_SUCCESS);
1213: }
1215: static PetscErrorCode MatDestroy_ScaLAPACK(Mat A)
1216: {
1217: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1218: Mat_ScaLAPACK_Grid *grid;
1219: PetscMPIInt iflg;
1220: MPI_Comm icomm;
1222: PetscFunctionBegin;
1223: PetscCall(MatStashDestroy_Private(&A->stash));
1224: PetscCall(PetscFree(a->loc));
1225: PetscCall(PetscFree(a->pivots));
1226: PetscCall(PetscCommDuplicate(PetscObjectComm((PetscObject)A), &icomm, NULL));
1227: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_ScaLAPACK_keyval, (void **)&grid, &iflg));
1228: if (--grid->grid_refct == 0) {
1229: Cblacs_freebuff(grid->ictxt, 1);
1230: Cblacs_gridexit(grid->ictxt);
1231: Cblacs_gridexit(grid->ictxrow);
1232: Cblacs_gridexit(grid->ictxcol);
1233: PetscCall(PetscFree(grid));
1234: PetscCallMPI(MPI_Comm_delete_attr(icomm, Petsc_ScaLAPACK_keyval));
1235: }
1236: PetscCall(PetscCommDestroy(&icomm));
1237: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", NULL));
1238: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
1239: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKSetBlockSizes_C", NULL));
1240: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKGetBlockSizes_C", NULL));
1241: PetscCall(PetscFree(A->data));
1242: PetscFunctionReturn(PETSC_SUCCESS);
1243: }
1245: static PetscErrorCode MatSetUp_ScaLAPACK(Mat A)
1246: {
1247: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1248: PetscBool flg;
1250: PetscFunctionBegin;
1251: PetscCall(PetscLayoutSetUp(A->rmap));
1252: PetscCall(PetscLayoutSetUp(A->cmap));
1254: /* check that the layout is as enforced by MatCreateScaLAPACK() */
1255: PetscCall(MatScaLAPACKCheckLayout(A->rmap, &flg));
1256: PetscCheck(flg, A->rmap->comm, PETSC_ERR_SUP, "MATSCALAPACK must have equal local row sizes in all processes (except possibly the last one), consider using MatCreateScaLAPACK");
1257: PetscCall(MatScaLAPACKCheckLayout(A->cmap, &flg));
1258: PetscCheck(flg, A->cmap->comm, PETSC_ERR_SUP, "MATSCALAPACK must have equal local column sizes in all processes (except possibly the last one), consider using MatCreateScaLAPACK");
1260: /* compute local sizes */
1261: PetscCall(PetscBLASIntCast(A->rmap->N, &a->M));
1262: PetscCall(PetscBLASIntCast(A->cmap->N, &a->N));
1263: a->locr = SCALAPACKnumroc_(&a->M, &a->mb, &a->grid->myrow, &a->rsrc, &a->grid->nprow);
1264: a->locc = SCALAPACKnumroc_(&a->N, &a->nb, &a->grid->mycol, &a->csrc, &a->grid->npcol);
1265: a->lld = PetscMax(1, a->locr);
1267: /* allocate local array */
1268: PetscCall(MatScaLAPACKSetPreallocation(A));
1270: /* set up ScaLAPACK descriptor */
1271: PetscCallScaLAPACKInfo("descinit", SCALAPACKdescinit_(a->desc, &a->M, &a->N, &a->mb, &a->nb, &a->rsrc, &a->csrc, &a->grid->ictxt, &a->lld, &info));
1272: PetscFunctionReturn(PETSC_SUCCESS);
1273: }
1275: static PetscErrorCode MatAssemblyBegin_ScaLAPACK(Mat A, MatAssemblyType type)
1276: {
1277: PetscInt nstash, reallocs;
1279: PetscFunctionBegin;
1280: if (A->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
1281: PetscCall(MatStashScatterBegin_Private(A, &A->stash, NULL));
1282: PetscCall(MatStashGetInfo_Private(&A->stash, &nstash, &reallocs));
1283: PetscCall(PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs));
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: static PetscErrorCode MatAssemblyEnd_ScaLAPACK(Mat A, MatAssemblyType type)
1288: {
1289: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1290: PetscMPIInt n;
1291: PetscInt i, flg, *row, *col;
1292: PetscScalar *val;
1293: PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc;
1295: PetscFunctionBegin;
1296: if (A->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
1297: while (1) {
1298: PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
1299: if (!flg) break;
1300: for (i = 0; i < n; i++) {
1301: PetscCall(PetscBLASIntCast(row[i] + 1, &gridx));
1302: PetscCall(PetscBLASIntCast(col[i] + 1, &gcidx));
1303: PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc));
1304: PetscCheck(rsrc == a->grid->myrow && csrc == a->grid->mycol, PetscObjectComm((PetscObject)A), PETSC_ERR_LIB, "Something went wrong, received value does not belong to this process");
1305: switch (A->insertmode) {
1306: case INSERT_VALUES:
1307: a->loc[lridx - 1 + (lcidx - 1) * a->lld] = val[i];
1308: break;
1309: case ADD_VALUES:
1310: a->loc[lridx - 1 + (lcidx - 1) * a->lld] += val[i];
1311: break;
1312: default:
1313: SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for InsertMode %d", (int)A->insertmode);
1314: }
1315: }
1316: }
1317: PetscCall(MatStashScatterEnd_Private(&A->stash));
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: static PetscErrorCode MatLoad_ScaLAPACK(Mat newMat, PetscViewer viewer)
1322: {
1323: Mat Adense, As;
1324: MPI_Comm comm;
1326: PetscFunctionBegin;
1327: PetscCall(PetscObjectGetComm((PetscObject)newMat, &comm));
1328: PetscCall(MatCreate(comm, &Adense));
1329: PetscCall(MatSetType(Adense, MATDENSE));
1330: PetscCall(MatLoad(Adense, viewer));
1331: PetscCall(MatConvert(Adense, MATSCALAPACK, MAT_INITIAL_MATRIX, &As));
1332: PetscCall(MatDestroy(&Adense));
1333: PetscCall(MatHeaderReplace(newMat, &As));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: static struct _MatOps MatOps_Values = {MatSetValues_ScaLAPACK,
1338: NULL,
1339: NULL,
1340: MatMult_ScaLAPACK,
1341: /* 4*/ MatMultAdd_ScaLAPACK,
1342: MatMultTranspose_ScaLAPACK,
1343: MatMultTransposeAdd_ScaLAPACK,
1344: MatSolve_ScaLAPACK,
1345: MatSolveAdd_ScaLAPACK,
1346: NULL,
1347: /*10*/ NULL,
1348: MatLUFactor_ScaLAPACK,
1349: MatCholeskyFactor_ScaLAPACK,
1350: NULL,
1351: MatTranspose_ScaLAPACK,
1352: /*15*/ MatGetInfo_ScaLAPACK,
1353: NULL,
1354: MatGetDiagonal_ScaLAPACK,
1355: MatDiagonalScale_ScaLAPACK,
1356: MatNorm_ScaLAPACK,
1357: /*20*/ MatAssemblyBegin_ScaLAPACK,
1358: MatAssemblyEnd_ScaLAPACK,
1359: MatSetOption_ScaLAPACK,
1360: MatZeroEntries_ScaLAPACK,
1361: /*24*/ NULL,
1362: MatLUFactorSymbolic_ScaLAPACK,
1363: MatLUFactorNumeric_ScaLAPACK,
1364: MatCholeskyFactorSymbolic_ScaLAPACK,
1365: MatCholeskyFactorNumeric_ScaLAPACK,
1366: /*29*/ MatSetUp_ScaLAPACK,
1367: NULL,
1368: NULL,
1369: NULL,
1370: NULL,
1371: /*34*/ MatDuplicate_ScaLAPACK,
1372: NULL,
1373: NULL,
1374: NULL,
1375: NULL,
1376: /*39*/ MatAXPY_ScaLAPACK,
1377: NULL,
1378: NULL,
1379: NULL,
1380: MatCopy_ScaLAPACK,
1381: /*44*/ NULL,
1382: MatScale_ScaLAPACK,
1383: MatShift_ScaLAPACK,
1384: NULL,
1385: NULL,
1386: /*49*/ NULL,
1387: NULL,
1388: NULL,
1389: NULL,
1390: NULL,
1391: /*54*/ NULL,
1392: NULL,
1393: NULL,
1394: NULL,
1395: NULL,
1396: /*59*/ NULL,
1397: MatDestroy_ScaLAPACK,
1398: MatView_ScaLAPACK,
1399: NULL,
1400: NULL,
1401: /*64*/ NULL,
1402: NULL,
1403: NULL,
1404: NULL,
1405: NULL,
1406: /*69*/ NULL,
1407: MatConvert_ScaLAPACK_Dense,
1408: NULL,
1409: NULL,
1410: NULL,
1411: /*74*/ NULL,
1412: NULL,
1413: NULL,
1414: NULL,
1415: MatLoad_ScaLAPACK,
1416: /*79*/ NULL,
1417: NULL,
1418: NULL,
1419: NULL,
1420: NULL,
1421: /*84*/ NULL,
1422: MatMatMultNumeric_ScaLAPACK,
1423: NULL,
1424: NULL,
1425: MatMatTransposeMultNumeric_ScaLAPACK,
1426: /*89*/ NULL,
1427: MatProductSetFromOptions_ScaLAPACK,
1428: NULL,
1429: NULL,
1430: MatConjugate_ScaLAPACK,
1431: /*94*/ NULL,
1432: NULL,
1433: NULL,
1434: NULL,
1435: NULL,
1436: /*99*/ NULL,
1437: MatMatSolve_ScaLAPACK,
1438: NULL,
1439: NULL,
1440: NULL,
1441: /*104*/ NULL,
1442: NULL,
1443: NULL,
1444: NULL,
1445: NULL,
1446: /*109*/ NULL,
1447: MatHermitianTranspose_ScaLAPACK,
1448: MatMultHermitianTranspose_ScaLAPACK,
1449: MatMultHermitianTransposeAdd_ScaLAPACK,
1450: NULL,
1451: /*114*/ NULL,
1452: NULL,
1453: NULL,
1454: NULL,
1455: NULL,
1456: /*119*/ NULL,
1457: MatTransposeMatMultNumeric_ScaLAPACK,
1458: NULL,
1459: NULL,
1460: /*124*/ NULL,
1461: NULL,
1462: NULL,
1463: NULL,
1464: NULL,
1465: /*129*/ NULL,
1466: NULL,
1467: NULL,
1468: NULL,
1469: NULL,
1470: /*134*/ NULL,
1471: NULL,
1472: NULL,
1473: NULL,
1474: NULL,
1475: NULL,
1476: /*140*/ NULL,
1477: NULL,
1478: NULL,
1479: NULL,
1480: NULL,
1481: /*144*/ NULL,
1482: NULL,
1483: NULL,
1484: NULL};
1486: static PetscErrorCode MatStashScatterBegin_ScaLAPACK(Mat mat, MatStash *stash, PetscInt *owners)
1487: {
1488: PetscInt *owner, *startv, *starti, bs2;
1489: PetscInt size = stash->size, nsends;
1490: PetscInt *sindices, **rindices, j, l;
1491: PetscScalar **rvalues, *svalues;
1492: MPI_Comm comm = stash->comm;
1493: MPI_Request *send_waits, *recv_waits, *recv_waits1, *recv_waits2;
1494: PetscMPIInt tag1 = stash->tag1, tag2 = stash->tag2, *sizes, *nlengths, nreceives, insends, ii;
1495: PetscInt *sp_idx, *sp_idy;
1496: PetscScalar *sp_val;
1497: PetscMatStashSpace space, space_next;
1498: PetscBLASInt gridx, gcidx, lridx, lcidx, rsrc, csrc;
1499: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)mat->data;
1501: PetscFunctionBegin;
1502: { /* make sure all processors are either in INSERTMODE or ADDMODE */
1503: InsertMode addv;
1504: PetscCallMPI(MPIU_Allreduce((PetscEnum *)&mat->insertmode, (PetscEnum *)&addv, 1, MPIU_ENUM, MPI_BOR, PetscObjectComm((PetscObject)mat)));
1505: PetscCheck(addv != (ADD_VALUES | INSERT_VALUES), PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Some processors inserted others added");
1506: mat->insertmode = addv; /* in case this processor had no cache */
1507: }
1509: bs2 = stash->bs * stash->bs;
1511: /* first count number of contributors to each processor */
1512: PetscCall(PetscCalloc1(size, &nlengths));
1513: PetscCall(PetscMalloc1(stash->n + 1, &owner));
1515: ii = j = 0;
1516: space = stash->space_head;
1517: while (space) {
1518: space_next = space->next;
1519: for (l = 0; l < space->local_used; l++) {
1520: PetscCall(PetscBLASIntCast(space->idx[l] + 1, &gridx));
1521: PetscCall(PetscBLASIntCast(space->idy[l] + 1, &gcidx));
1522: PetscCallBLAS("SCALAPACKinfog2l", SCALAPACKinfog2l_(&gridx, &gcidx, a->desc, &a->grid->nprow, &a->grid->npcol, &a->grid->myrow, &a->grid->mycol, &lridx, &lcidx, &rsrc, &csrc));
1523: j = Cblacs_pnum(a->grid->ictxt, rsrc, csrc);
1524: nlengths[j]++;
1525: owner[ii] = j;
1526: ii++;
1527: }
1528: space = space_next;
1529: }
1531: /* Now check what procs get messages - and compute nsends. */
1532: PetscCall(PetscCalloc1(size, &sizes));
1533: nsends = 0;
1534: for (PetscMPIInt i = 0; i < size; i++) {
1535: if (nlengths[i]) {
1536: sizes[i] = 1;
1537: nsends++;
1538: }
1539: }
1541: {
1542: PetscMPIInt *onodes, *olengths;
1544: /* Determine the number of messages to expect, their lengths, from from-ids */
1545: PetscCall(PetscGatherNumberOfMessages(comm, sizes, nlengths, &nreceives));
1546: PetscCall(PetscMPIIntCast(nsends, &insends));
1547: PetscCall(PetscGatherMessageLengths(comm, insends, nreceives, nlengths, &onodes, &olengths));
1548: /* since clubbing row,col - lengths are multiplied by 2 */
1549: for (PetscMPIInt i = 0; i < nreceives; i++) olengths[i] *= 2;
1550: PetscCall(PetscPostIrecvInt(comm, tag1, nreceives, onodes, olengths, &rindices, &recv_waits1));
1551: /* values are size 'bs2' lengths (and remove earlier factor 2 */
1552: for (PetscMPIInt i = 0; i < nreceives; i++) olengths[i] = (PetscMPIInt)(olengths[i] * bs2 / 2);
1553: PetscCall(PetscPostIrecvScalar(comm, tag2, nreceives, onodes, olengths, &rvalues, &recv_waits2));
1554: PetscCall(PetscFree(onodes));
1555: PetscCall(PetscFree(olengths));
1556: }
1558: /* do sends:
1559: 1) starts[i] gives the starting index in svalues for stuff going to
1560: the ith processor
1561: */
1562: PetscCall(PetscMalloc2(bs2 * stash->n, &svalues, 2 * (stash->n + 1), &sindices));
1563: PetscCall(PetscMalloc1(2 * nsends, &send_waits));
1564: PetscCall(PetscMalloc2(size, &startv, size, &starti));
1565: /* use 2 sends the first with all_a, the next with all_i and all_j */
1566: startv[0] = 0;
1567: starti[0] = 0;
1568: for (PetscMPIInt i = 1; i < size; i++) {
1569: startv[i] = startv[i - 1] + nlengths[i - 1];
1570: starti[i] = starti[i - 1] + 2 * nlengths[i - 1];
1571: }
1573: ii = 0;
1574: space = stash->space_head;
1575: while (space) {
1576: space_next = space->next;
1577: sp_idx = space->idx;
1578: sp_idy = space->idy;
1579: sp_val = space->val;
1580: for (l = 0; l < space->local_used; l++) {
1581: j = owner[ii];
1582: if (bs2 == 1) {
1583: svalues[startv[j]] = sp_val[l];
1584: } else {
1585: PetscScalar *buf1, *buf2;
1586: buf1 = svalues + bs2 * startv[j];
1587: buf2 = space->val + bs2 * l;
1588: for (PetscInt k = 0; k < bs2; k++) buf1[k] = buf2[k];
1589: }
1590: sindices[starti[j]] = sp_idx[l];
1591: sindices[starti[j] + nlengths[j]] = sp_idy[l];
1592: startv[j]++;
1593: starti[j]++;
1594: ii++;
1595: }
1596: space = space_next;
1597: }
1598: startv[0] = 0;
1599: for (PetscMPIInt i = 1; i < size; i++) startv[i] = startv[i - 1] + nlengths[i - 1];
1601: for (PetscMPIInt i = 0, count = 0; i < size; i++) {
1602: if (sizes[i]) {
1603: PetscCallMPI(MPIU_Isend(sindices + 2 * startv[i], 2 * nlengths[i], MPIU_INT, i, tag1, comm, send_waits + count++));
1604: PetscCallMPI(MPIU_Isend(svalues + bs2 * startv[i], bs2 * nlengths[i], MPIU_SCALAR, i, tag2, comm, send_waits + count++));
1605: }
1606: }
1607: if (PetscDefined(USE_INFO)) {
1608: PetscCall(PetscInfo(NULL, "No of messages: %" PetscInt_FMT "\n", nsends));
1609: for (PetscMPIInt i = 0; i < size; i++) {
1610: if (sizes[i]) PetscCall(PetscInfo(NULL, "Mesg_to: %d: size: %zu bytes\n", i, (size_t)(nlengths[i] * (bs2 * sizeof(PetscScalar) + 2 * sizeof(PetscInt)))));
1611: }
1612: }
1613: PetscCall(PetscFree(nlengths));
1614: PetscCall(PetscFree(owner));
1615: PetscCall(PetscFree2(startv, starti));
1616: PetscCall(PetscFree(sizes));
1618: /* recv_waits need to be contiguous for MatStashScatterGetMesg_Private() */
1619: PetscCall(PetscMalloc1(2 * nreceives, &recv_waits));
1621: for (PetscMPIInt i = 0; i < nreceives; i++) {
1622: recv_waits[2 * i] = recv_waits1[i];
1623: recv_waits[2 * i + 1] = recv_waits2[i];
1624: }
1625: stash->recv_waits = recv_waits;
1627: PetscCall(PetscFree(recv_waits1));
1628: PetscCall(PetscFree(recv_waits2));
1630: stash->svalues = svalues;
1631: stash->sindices = sindices;
1632: stash->rvalues = rvalues;
1633: stash->rindices = rindices;
1634: stash->send_waits = send_waits;
1635: stash->nsends = (PetscMPIInt)nsends;
1636: stash->nrecvs = nreceives;
1637: stash->reproduce_count = 0;
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: static PetscErrorCode MatScaLAPACKSetBlockSizes_ScaLAPACK(Mat A, PetscInt mb, PetscInt nb)
1642: {
1643: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1645: PetscFunctionBegin;
1646: PetscCheck(!A->preallocated, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Cannot change block sizes after MatSetUp");
1647: PetscCheck(mb >= 1 || mb == PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "mb %" PetscInt_FMT " must be at least 1", mb);
1648: PetscCheck(nb >= 1 || nb == PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "nb %" PetscInt_FMT " must be at least 1", nb);
1649: PetscCall(PetscBLASIntCast((mb == PETSC_DECIDE) ? DEFAULT_BLOCKSIZE : mb, &a->mb));
1650: PetscCall(PetscBLASIntCast((nb == PETSC_DECIDE) ? a->mb : nb, &a->nb));
1651: PetscFunctionReturn(PETSC_SUCCESS);
1652: }
1654: /*@
1655: MatScaLAPACKSetBlockSizes - Sets the block sizes to be used for the distribution of
1656: the `MATSCALAPACK` matrix
1658: Logically Collective
1660: Input Parameters:
1661: + A - a `MATSCALAPACK` matrix
1662: . mb - the row block size
1663: - nb - the column block size
1665: Level: intermediate
1667: Note:
1668: This block size has a different meaning from the block size associated with `MatSetBlockSize()` used for sparse matrices
1670: .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MatCreateScaLAPACK()`, `MatScaLAPACKGetBlockSizes()`
1671: @*/
1672: PetscErrorCode MatScaLAPACKSetBlockSizes(Mat A, PetscInt mb, PetscInt nb)
1673: {
1674: PetscFunctionBegin;
1678: PetscTryMethod(A, "MatScaLAPACKSetBlockSizes_C", (Mat, PetscInt, PetscInt), (A, mb, nb));
1679: PetscFunctionReturn(PETSC_SUCCESS);
1680: }
1682: static PetscErrorCode MatScaLAPACKGetBlockSizes_ScaLAPACK(Mat A, PetscInt *mb, PetscInt *nb)
1683: {
1684: Mat_ScaLAPACK *a = (Mat_ScaLAPACK *)A->data;
1686: PetscFunctionBegin;
1687: if (mb) *mb = a->mb;
1688: if (nb) *nb = a->nb;
1689: PetscFunctionReturn(PETSC_SUCCESS);
1690: }
1692: /*@
1693: MatScaLAPACKGetBlockSizes - Gets the block sizes used in the distribution of
1694: the `MATSCALAPACK` matrix
1696: Not Collective
1698: Input Parameter:
1699: . A - a `MATSCALAPACK` matrix
1701: Output Parameters:
1702: + mb - the row block size
1703: - nb - the column block size
1705: Level: intermediate
1707: Note:
1708: This block size has a different meaning from the block size associated with `MatSetBlockSize()` used for sparse matrices
1710: .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MatCreateScaLAPACK()`, `MatScaLAPACKSetBlockSizes()`
1711: @*/
1712: PetscErrorCode MatScaLAPACKGetBlockSizes(Mat A, PetscInt *mb, PetscInt *nb)
1713: {
1714: PetscFunctionBegin;
1716: PetscUseMethod(A, "MatScaLAPACKGetBlockSizes_C", (Mat, PetscInt *, PetscInt *), (A, mb, nb));
1717: PetscFunctionReturn(PETSC_SUCCESS);
1718: }
1720: PETSC_INTERN PetscErrorCode MatStashScatterGetMesg_Ref(MatStash *, PetscMPIInt *, PetscInt **, PetscInt **, PetscScalar **, PetscInt *);
1721: PETSC_INTERN PetscErrorCode MatStashScatterEnd_Ref(MatStash *);
1723: /*MC
1724: MATSCALAPACK = "scalapack" - A matrix type for dense matrices using the ScaLAPACK package
1726: Use `./configure --download-scalapack` to install PETSc to use ScaLAPACK
1728: Options Database Keys:
1729: + -mat_type scalapack - sets the matrix type to `MATSCALAPACK`
1730: . -pc_factor_mat_solver_type scalapack - to use this direct solver with the option `-pc_type lu`
1731: . -mat_scalapack_grid_height - sets Grid Height for 2D cyclic ordering of internal matrix
1732: - -mat_scalapack_block_sizes - size of the blocks to use (one or two integers separated by comma)
1734: Level: intermediate
1736: Note:
1737: Note unlike most matrix formats, this format does not store all the matrix entries for a contiguous
1738: range of rows on an MPI rank. Use `MatGetOwnershipIS()` to determine what values are stored on
1739: the given rank.
1741: .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MATDENSE`, `MATELEMENTAL`, `MatGetOwnershipIS()`, `MatCreateScaLAPACK()`
1742: M*/
1744: PETSC_EXTERN PetscErrorCode MatCreate_ScaLAPACK(Mat A)
1745: {
1746: Mat_ScaLAPACK *a;
1747: PetscBool flg;
1748: PetscMPIInt iflg;
1749: Mat_ScaLAPACK_Grid *grid;
1750: MPI_Comm icomm;
1751: PetscBLASInt nprow, npcol, myrow, mycol;
1752: PetscInt optv1, k = 2, array[2] = {0, 0};
1753: PetscMPIInt size;
1755: PetscFunctionBegin;
1756: A->ops[0] = MatOps_Values;
1757: A->insertmode = NOT_SET_VALUES;
1759: PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)A), 1, &A->stash));
1760: A->stash.ScatterBegin = MatStashScatterBegin_ScaLAPACK;
1761: A->stash.ScatterGetMesg = MatStashScatterGetMesg_Ref;
1762: A->stash.ScatterEnd = MatStashScatterEnd_Ref;
1763: A->stash.ScatterDestroy = NULL;
1765: PetscCall(PetscNew(&a));
1766: A->data = (void *)a;
1768: /* Grid needs to be shared between multiple Mats on the same communicator, implement by attribute caching on the MPI_Comm */
1769: if (Petsc_ScaLAPACK_keyval == MPI_KEYVAL_INVALID) {
1770: PetscCallMPI(MPI_Comm_create_keyval(MPI_COMM_NULL_COPY_FN, MPI_COMM_NULL_DELETE_FN, &Petsc_ScaLAPACK_keyval, NULL));
1771: PetscCall(PetscRegisterFinalize(Petsc_ScaLAPACK_keyval_free));
1772: PetscCall(PetscCitationsRegister(ScaLAPACKCitation, &ScaLAPACKCite));
1773: }
1774: PetscCall(PetscCommDuplicate(PetscObjectComm((PetscObject)A), &icomm, NULL));
1775: PetscCallMPI(MPI_Comm_get_attr(icomm, Petsc_ScaLAPACK_keyval, (void **)&grid, &iflg));
1776: if (!iflg) {
1777: PetscCall(PetscNew(&grid));
1779: PetscCallMPI(MPI_Comm_size(icomm, &size));
1780: PetscCall(PetscBLASIntCast(PetscSqrtReal((PetscReal)size) + 0.001, &grid->nprow));
1782: PetscOptionsBegin(PetscObjectComm((PetscObject)A), ((PetscObject)A)->prefix, "ScaLAPACK Grid Options", "Mat");
1783: PetscCall(PetscOptionsInt("-mat_scalapack_grid_height", "Grid Height", "None", grid->nprow, &optv1, &flg));
1784: if (flg) {
1785: PetscCheck(size % optv1 == 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Grid Height %" PetscInt_FMT " must evenly divide CommSize %d", optv1, size);
1786: PetscCall(PetscBLASIntCast(optv1, &grid->nprow));
1787: }
1788: PetscOptionsEnd();
1790: if (size % grid->nprow) grid->nprow = 1; /* cannot use a squarish grid, use a 1d grid */
1791: grid->npcol = size / grid->nprow;
1792: PetscCall(PetscBLASIntCast(grid->nprow, &nprow));
1793: PetscCall(PetscBLASIntCast(grid->npcol, &npcol));
1794: grid->ictxt = Csys2blacs_handle(icomm);
1795: Cblacs_gridinit(&grid->ictxt, "R", nprow, npcol);
1796: Cblacs_gridinfo(grid->ictxt, &nprow, &npcol, &myrow, &mycol);
1797: grid->grid_refct = 1;
1798: grid->nprow = nprow;
1799: grid->npcol = npcol;
1800: grid->myrow = myrow;
1801: grid->mycol = mycol;
1802: /* auxiliary 1d BLACS contexts for 1xsize and sizex1 grids */
1803: grid->ictxrow = Csys2blacs_handle(icomm);
1804: Cblacs_gridinit(&grid->ictxrow, "R", 1, size);
1805: grid->ictxcol = Csys2blacs_handle(icomm);
1806: Cblacs_gridinit(&grid->ictxcol, "R", size, 1);
1807: PetscCallMPI(MPI_Comm_set_attr(icomm, Petsc_ScaLAPACK_keyval, (void *)grid));
1809: } else grid->grid_refct++;
1810: PetscCall(PetscCommDestroy(&icomm));
1811: a->grid = grid;
1812: a->mb = DEFAULT_BLOCKSIZE;
1813: a->nb = DEFAULT_BLOCKSIZE;
1815: PetscOptionsBegin(PetscObjectComm((PetscObject)A), NULL, "ScaLAPACK Options", "Mat");
1816: PetscCall(PetscOptionsIntArray("-mat_scalapack_block_sizes", "Size of the blocks to use (one or two comma-separated integers)", "MatCreateScaLAPACK", array, &k, &flg));
1817: if (flg) {
1818: a->mb = (PetscMPIInt)array[0];
1819: a->nb = (k > 1) ? (PetscMPIInt)array[1] : a->mb;
1820: }
1821: PetscOptionsEnd();
1823: a->roworiented = PETSC_TRUE;
1824: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatGetOwnershipIS_C", MatGetOwnershipIS_ScaLAPACK));
1825: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKSetBlockSizes_C", MatScaLAPACKSetBlockSizes_ScaLAPACK));
1826: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatScaLAPACKGetBlockSizes_C", MatScaLAPACKGetBlockSizes_ScaLAPACK));
1827: PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSCALAPACK));
1828: PetscFunctionReturn(PETSC_SUCCESS);
1829: }
1831: /*@
1832: MatCreateScaLAPACK - Creates a dense parallel matrix in ScaLAPACK format
1833: (2D block cyclic distribution) for a `MATSCALAPACK` matrix
1835: Collective
1837: Input Parameters:
1838: + comm - MPI communicator
1839: . mb - row block size (or `PETSC_DECIDE` to have it set)
1840: . nb - column block size (or `PETSC_DECIDE` to have it set)
1841: . M - number of global rows
1842: . N - number of global columns
1843: . rsrc - coordinate of process that owns the first row of the distributed matrix
1844: - csrc - coordinate of process that owns the first column of the distributed matrix
1846: Output Parameter:
1847: . A - the matrix
1849: Options Database Key:
1850: . -mat_scalapack_block_sizes - size of the blocks to use (one or two integers separated by comma)
1852: Level: intermediate
1854: Notes:
1855: If `PETSC_DECIDE` is used for the block sizes, then an appropriate value is chosen
1857: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
1858: MatXXXXSetPreallocation() paradigm instead of this routine directly.
1859: [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]
1861: Storage is completely managed by ScaLAPACK, so this requires PETSc to be
1862: configured with ScaLAPACK. In particular, PETSc's local sizes lose
1863: significance and are thus ignored. The block sizes refer to the values
1864: used for the distributed matrix, not the same meaning as in `MATBAIJ`.
1866: .seealso: [](ch_matrices), `Mat`, `MATSCALAPACK`, `MATDENSE`, `MATELEMENTAL`, `MatCreate()`, `MatCreateDense()`, `MatSetValues()`
1867: @*/
1868: PetscErrorCode MatCreateScaLAPACK(MPI_Comm comm, PetscInt mb, PetscInt nb, PetscInt M, PetscInt N, PetscInt rsrc, PetscInt csrc, Mat *A)
1869: {
1870: Mat_ScaLAPACK *a;
1871: PetscInt m, n;
1873: PetscFunctionBegin;
1874: PetscCall(MatCreate(comm, A));
1875: PetscCall(MatSetType(*A, MATSCALAPACK));
1876: PetscCheck(M != PETSC_DECIDE && N != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot use PETSC_DECIDE for matrix dimensions");
1877: /* rows and columns are NOT distributed according to PetscSplitOwnership */
1878: m = PETSC_DECIDE;
1879: PetscCall(PetscSplitOwnershipEqual(comm, &m, &M));
1880: n = PETSC_DECIDE;
1881: PetscCall(PetscSplitOwnershipEqual(comm, &n, &N));
1882: PetscCall(MatSetSizes(*A, m, n, M, N));
1883: a = (Mat_ScaLAPACK *)(*A)->data;
1884: PetscCall(PetscBLASIntCast(M, &a->M));
1885: PetscCall(PetscBLASIntCast(N, &a->N));
1886: PetscCall(PetscBLASIntCast((mb == PETSC_DECIDE) ? DEFAULT_BLOCKSIZE : mb, &a->mb));
1887: PetscCall(PetscBLASIntCast((nb == PETSC_DECIDE) ? a->mb : nb, &a->nb));
1888: PetscCall(PetscBLASIntCast(rsrc, &a->rsrc));
1889: PetscCall(PetscBLASIntCast(csrc, &a->csrc));
1890: PetscCall(MatSetUp(*A));
1891: PetscFunctionReturn(PETSC_SUCCESS);
1892: }