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