1 2 /* 3 Basic functions for basic parallel dense matrices. 4 */ 5 6 7 #include <../src/mat/impls/dense/mpi/mpidense.h> /*I "petscmat.h" I*/ 8 #if defined(PETSC_HAVE_PLAPACK) 9 static PetscMPIInt Plapack_nprows,Plapack_npcols,Plapack_ierror,Plapack_nb_alg; 10 static MPI_Comm Plapack_comm_2d; 11 EXTERN_C_BEGIN 12 #include <PLA.h> 13 EXTERN_C_END 14 15 typedef struct { 16 PLA_Obj A,pivots; 17 PLA_Template templ; 18 MPI_Datatype datatype; 19 PetscInt nb,rstart; 20 VecScatter ctx; 21 IS is_pla,is_petsc; 22 PetscBool pla_solved; 23 MatStructure mstruct; 24 } Mat_Plapack; 25 #endif 26 27 #undef __FUNCT__ 28 #define __FUNCT__ "MatDenseGetLocalMatrix" 29 /*@ 30 31 MatDenseGetLocalMatrix - For a MATMPIDENSE or MATSEQDENSE matrix returns the sequential 32 matrix that represents the operator. For sequential matrices it returns itself. 33 34 Input Parameter: 35 . A - the Seq or MPI dense matrix 36 37 Output Parameter: 38 . B - the inner matrix 39 40 Level: intermediate 41 42 @*/ 43 PetscErrorCode MatDenseGetLocalMatrix(Mat A,Mat *B) 44 { 45 Mat_MPIDense *mat = (Mat_MPIDense*)A->data; 46 PetscErrorCode ierr; 47 PetscBool flg; 48 49 PetscFunctionBegin; 50 ierr = PetscTypeCompare((PetscObject)A,MATMPIDENSE,&flg);CHKERRQ(ierr); 51 if (flg) { 52 *B = mat->A; 53 } else { 54 *B = A; 55 } 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "MatGetRow_MPIDense" 61 PetscErrorCode MatGetRow_MPIDense(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 62 { 63 Mat_MPIDense *mat = (Mat_MPIDense*)A->data; 64 PetscErrorCode ierr; 65 PetscInt lrow,rstart = A->rmap->rstart,rend = A->rmap->rend; 66 67 PetscFunctionBegin; 68 if (row < rstart || row >= rend) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"only local rows"); 69 lrow = row - rstart; 70 ierr = MatGetRow(mat->A,lrow,nz,(const PetscInt **)idx,(const PetscScalar **)v);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 #undef __FUNCT__ 75 #define __FUNCT__ "MatRestoreRow_MPIDense" 76 PetscErrorCode MatRestoreRow_MPIDense(Mat mat,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 77 { 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 82 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 83 PetscFunctionReturn(0); 84 } 85 86 EXTERN_C_BEGIN 87 #undef __FUNCT__ 88 #define __FUNCT__ "MatGetDiagonalBlock_MPIDense" 89 PetscErrorCode MatGetDiagonalBlock_MPIDense(Mat A,Mat *a) 90 { 91 Mat_MPIDense *mdn = (Mat_MPIDense*)A->data; 92 PetscErrorCode ierr; 93 PetscInt m = A->rmap->n,rstart = A->rmap->rstart; 94 PetscScalar *array; 95 MPI_Comm comm; 96 Mat B; 97 98 PetscFunctionBegin; 99 if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only square matrices supported."); 100 101 ierr = PetscObjectQuery((PetscObject)A,"DiagonalBlock",(PetscObject*)&B);CHKERRQ(ierr); 102 if (!B) { 103 ierr = PetscObjectGetComm((PetscObject)(mdn->A),&comm);CHKERRQ(ierr); 104 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 105 ierr = MatSetSizes(B,m,m,m,m);CHKERRQ(ierr); 106 ierr = MatSetType(B,((PetscObject)mdn->A)->type_name);CHKERRQ(ierr); 107 ierr = MatGetArray(mdn->A,&array);CHKERRQ(ierr); 108 ierr = MatSeqDenseSetPreallocation(B,array+m*rstart);CHKERRQ(ierr); 109 ierr = MatRestoreArray(mdn->A,&array);CHKERRQ(ierr); 110 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 111 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 112 ierr = PetscObjectCompose((PetscObject)A,"DiagonalBlock",(PetscObject)B);CHKERRQ(ierr); 113 *a = B; 114 ierr = MatDestroy(&B);CHKERRQ(ierr); 115 } else { 116 *a = B; 117 } 118 PetscFunctionReturn(0); 119 } 120 EXTERN_C_END 121 122 #undef __FUNCT__ 123 #define __FUNCT__ "MatSetValues_MPIDense" 124 PetscErrorCode MatSetValues_MPIDense(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],const PetscScalar v[],InsertMode addv) 125 { 126 Mat_MPIDense *A = (Mat_MPIDense*)mat->data; 127 PetscErrorCode ierr; 128 PetscInt i,j,rstart = mat->rmap->rstart,rend = mat->rmap->rend,row; 129 PetscBool roworiented = A->roworiented; 130 131 PetscFunctionBegin; 132 if (v) PetscValidScalarPointer(v,6); 133 for (i=0; i<m; i++) { 134 if (idxm[i] < 0) continue; 135 if (idxm[i] >= mat->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 136 if (idxm[i] >= rstart && idxm[i] < rend) { 137 row = idxm[i] - rstart; 138 if (roworiented) { 139 ierr = MatSetValues(A->A,1,&row,n,idxn,v+i*n,addv);CHKERRQ(ierr); 140 } else { 141 for (j=0; j<n; j++) { 142 if (idxn[j] < 0) continue; 143 if (idxn[j] >= mat->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 144 ierr = MatSetValues(A->A,1,&row,1,&idxn[j],v+i+j*m,addv);CHKERRQ(ierr); 145 } 146 } 147 } else { 148 if (!A->donotstash) { 149 if (roworiented) { 150 ierr = MatStashValuesRow_Private(&mat->stash,idxm[i],n,idxn,v+i*n,PETSC_FALSE);CHKERRQ(ierr); 151 } else { 152 ierr = MatStashValuesCol_Private(&mat->stash,idxm[i],n,idxn,v+i,m,PETSC_FALSE);CHKERRQ(ierr); 153 } 154 } 155 } 156 } 157 PetscFunctionReturn(0); 158 } 159 160 #undef __FUNCT__ 161 #define __FUNCT__ "MatGetValues_MPIDense" 162 PetscErrorCode MatGetValues_MPIDense(Mat mat,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[]) 163 { 164 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 165 PetscErrorCode ierr; 166 PetscInt i,j,rstart = mat->rmap->rstart,rend = mat->rmap->rend,row; 167 168 PetscFunctionBegin; 169 for (i=0; i<m; i++) { 170 if (idxm[i] < 0) continue; /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row"); */ 171 if (idxm[i] >= mat->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 172 if (idxm[i] >= rstart && idxm[i] < rend) { 173 row = idxm[i] - rstart; 174 for (j=0; j<n; j++) { 175 if (idxn[j] < 0) continue; /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column"); */ 176 if (idxn[j] >= mat->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 177 ierr = MatGetValues(mdn->A,1,&row,1,&idxn[j],v+i*n+j);CHKERRQ(ierr); 178 } 179 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only local values currently supported"); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "MatGetArray_MPIDense" 186 PetscErrorCode MatGetArray_MPIDense(Mat A,PetscScalar *array[]) 187 { 188 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 ierr = MatGetArray(a->A,array);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "MatGetSubMatrix_MPIDense" 198 static PetscErrorCode MatGetSubMatrix_MPIDense(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B) 199 { 200 Mat_MPIDense *mat = (Mat_MPIDense*)A->data,*newmatd; 201 Mat_SeqDense *lmat = (Mat_SeqDense*)mat->A->data; 202 PetscErrorCode ierr; 203 PetscInt i,j,rstart,rend,nrows,ncols,Ncols,nlrows,nlcols; 204 const PetscInt *irow,*icol; 205 PetscScalar *av,*bv,*v = lmat->v; 206 Mat newmat; 207 IS iscol_local; 208 209 PetscFunctionBegin; 210 ierr = ISAllGather(iscol,&iscol_local);CHKERRQ(ierr); 211 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 212 ierr = ISGetIndices(iscol_local,&icol);CHKERRQ(ierr); 213 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 214 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 215 ierr = ISGetSize(iscol,&Ncols);CHKERRQ(ierr); /* global number of columns, size of iscol_local */ 216 217 /* No parallel redistribution currently supported! Should really check each index set 218 to comfirm that it is OK. ... Currently supports only submatrix same partitioning as 219 original matrix! */ 220 221 ierr = MatGetLocalSize(A,&nlrows,&nlcols);CHKERRQ(ierr); 222 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 223 224 /* Check submatrix call */ 225 if (scall == MAT_REUSE_MATRIX) { 226 /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); */ 227 /* Really need to test rows and column sizes! */ 228 newmat = *B; 229 } else { 230 /* Create and fill new matrix */ 231 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 232 ierr = MatSetSizes(newmat,nrows,ncols,PETSC_DECIDE,Ncols);CHKERRQ(ierr); 233 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 234 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 235 } 236 237 /* Now extract the data pointers and do the copy, column at a time */ 238 newmatd = (Mat_MPIDense*)newmat->data; 239 bv = ((Mat_SeqDense *)newmatd->A->data)->v; 240 241 for (i=0; i<Ncols; i++) { 242 av = v + ((Mat_SeqDense *)mat->A->data)->lda*icol[i]; 243 for (j=0; j<nrows; j++) { 244 *bv++ = av[irow[j] - rstart]; 245 } 246 } 247 248 /* Assemble the matrices so that the correct flags are set */ 249 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 250 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 251 252 /* Free work space */ 253 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 254 ierr = ISRestoreIndices(iscol_local,&icol);CHKERRQ(ierr); 255 *B = newmat; 256 PetscFunctionReturn(0); 257 } 258 259 #undef __FUNCT__ 260 #define __FUNCT__ "MatRestoreArray_MPIDense" 261 PetscErrorCode MatRestoreArray_MPIDense(Mat A,PetscScalar *array[]) 262 { 263 PetscFunctionBegin; 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "MatAssemblyBegin_MPIDense" 269 PetscErrorCode MatAssemblyBegin_MPIDense(Mat mat,MatAssemblyType mode) 270 { 271 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 272 MPI_Comm comm = ((PetscObject)mat)->comm; 273 PetscErrorCode ierr; 274 PetscInt nstash,reallocs; 275 InsertMode addv; 276 277 PetscFunctionBegin; 278 /* make sure all processors are either in INSERTMODE or ADDMODE */ 279 ierr = MPI_Allreduce(&mat->insertmode,&addv,1,MPI_INT,MPI_BOR,comm);CHKERRQ(ierr); 280 if (addv == (ADD_VALUES|INSERT_VALUES)) SETERRQ(((PetscObject)mat)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot mix adds/inserts on different procs"); 281 mat->insertmode = addv; /* in case this processor had no cache */ 282 283 ierr = MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);CHKERRQ(ierr); 284 ierr = MatStashGetInfo_Private(&mat->stash,&nstash,&reallocs);CHKERRQ(ierr); 285 ierr = PetscInfo2(mdn->A,"Stash has %D entries, uses %D mallocs.\n",nstash,reallocs);CHKERRQ(ierr); 286 PetscFunctionReturn(0); 287 } 288 289 #undef __FUNCT__ 290 #define __FUNCT__ "MatAssemblyEnd_MPIDense" 291 PetscErrorCode MatAssemblyEnd_MPIDense(Mat mat,MatAssemblyType mode) 292 { 293 Mat_MPIDense *mdn=(Mat_MPIDense*)mat->data; 294 PetscErrorCode ierr; 295 PetscInt i,*row,*col,flg,j,rstart,ncols; 296 PetscMPIInt n; 297 PetscScalar *val; 298 InsertMode addv=mat->insertmode; 299 300 PetscFunctionBegin; 301 /* wait on receives */ 302 while (1) { 303 ierr = MatStashScatterGetMesg_Private(&mat->stash,&n,&row,&col,&val,&flg);CHKERRQ(ierr); 304 if (!flg) break; 305 306 for (i=0; i<n;) { 307 /* Now identify the consecutive vals belonging to the same row */ 308 for (j=i,rstart=row[j]; j<n; j++) { if (row[j] != rstart) break; } 309 if (j < n) ncols = j-i; 310 else ncols = n-i; 311 /* Now assemble all these values with a single function call */ 312 ierr = MatSetValues_MPIDense(mat,1,row+i,ncols,col+i,val+i,addv);CHKERRQ(ierr); 313 i = j; 314 } 315 } 316 ierr = MatStashScatterEnd_Private(&mat->stash);CHKERRQ(ierr); 317 318 ierr = MatAssemblyBegin(mdn->A,mode);CHKERRQ(ierr); 319 ierr = MatAssemblyEnd(mdn->A,mode);CHKERRQ(ierr); 320 321 if (!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) { 322 ierr = MatSetUpMultiply_MPIDense(mat);CHKERRQ(ierr); 323 } 324 PetscFunctionReturn(0); 325 } 326 327 #undef __FUNCT__ 328 #define __FUNCT__ "MatZeroEntries_MPIDense" 329 PetscErrorCode MatZeroEntries_MPIDense(Mat A) 330 { 331 PetscErrorCode ierr; 332 Mat_MPIDense *l = (Mat_MPIDense*)A->data; 333 334 PetscFunctionBegin; 335 ierr = MatZeroEntries(l->A);CHKERRQ(ierr); 336 PetscFunctionReturn(0); 337 } 338 339 /* the code does not do the diagonal entries correctly unless the 340 matrix is square and the column and row owerships are identical. 341 This is a BUG. The only way to fix it seems to be to access 342 mdn->A and mdn->B directly and not through the MatZeroRows() 343 routine. 344 */ 345 #undef __FUNCT__ 346 #define __FUNCT__ "MatZeroRows_MPIDense" 347 PetscErrorCode MatZeroRows_MPIDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 348 { 349 Mat_MPIDense *l = (Mat_MPIDense*)A->data; 350 PetscErrorCode ierr; 351 PetscInt i,*owners = A->rmap->range; 352 PetscInt *nprocs,j,idx,nsends; 353 PetscInt nmax,*svalues,*starts,*owner,nrecvs; 354 PetscInt *rvalues,tag = ((PetscObject)A)->tag,count,base,slen,*source; 355 PetscInt *lens,*lrows,*values; 356 PetscMPIInt n,imdex,rank = l->rank,size = l->size; 357 MPI_Comm comm = ((PetscObject)A)->comm; 358 MPI_Request *send_waits,*recv_waits; 359 MPI_Status recv_status,*send_status; 360 PetscBool found; 361 const PetscScalar *xx; 362 PetscScalar *bb; 363 364 PetscFunctionBegin; 365 if (A->rmap->N != A->cmap->N) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Only handles square matrices"); 366 if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only handles matrices with identical column and row ownership"); 367 /* first count number of contributors to each processor */ 368 ierr = PetscMalloc(2*size*sizeof(PetscInt),&nprocs);CHKERRQ(ierr); 369 ierr = PetscMemzero(nprocs,2*size*sizeof(PetscInt));CHKERRQ(ierr); 370 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&owner);CHKERRQ(ierr); /* see note*/ 371 for (i=0; i<N; i++) { 372 idx = rows[i]; 373 found = PETSC_FALSE; 374 for (j=0; j<size; j++) { 375 if (idx >= owners[j] && idx < owners[j+1]) { 376 nprocs[2*j]++; nprocs[2*j+1] = 1; owner[i] = j; found = PETSC_TRUE; break; 377 } 378 } 379 if (!found) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Index out of range"); 380 } 381 nsends = 0; for (i=0; i<size; i++) { nsends += nprocs[2*i+1];} 382 383 /* inform other processors of number of messages and max length*/ 384 ierr = PetscMaxSum(comm,nprocs,&nmax,&nrecvs);CHKERRQ(ierr); 385 386 /* post receives: */ 387 ierr = PetscMalloc((nrecvs+1)*(nmax+1)*sizeof(PetscInt),&rvalues);CHKERRQ(ierr); 388 ierr = PetscMalloc((nrecvs+1)*sizeof(MPI_Request),&recv_waits);CHKERRQ(ierr); 389 for (i=0; i<nrecvs; i++) { 390 ierr = MPI_Irecv(rvalues+nmax*i,nmax,MPIU_INT,MPI_ANY_SOURCE,tag,comm,recv_waits+i);CHKERRQ(ierr); 391 } 392 393 /* do sends: 394 1) starts[i] gives the starting index in svalues for stuff going to 395 the ith processor 396 */ 397 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&svalues);CHKERRQ(ierr); 398 ierr = PetscMalloc((nsends+1)*sizeof(MPI_Request),&send_waits);CHKERRQ(ierr); 399 ierr = PetscMalloc((size+1)*sizeof(PetscInt),&starts);CHKERRQ(ierr); 400 starts[0] = 0; 401 for (i=1; i<size; i++) { starts[i] = starts[i-1] + nprocs[2*i-2];} 402 for (i=0; i<N; i++) { 403 svalues[starts[owner[i]]++] = rows[i]; 404 } 405 406 starts[0] = 0; 407 for (i=1; i<size+1; i++) { starts[i] = starts[i-1] + nprocs[2*i-2];} 408 count = 0; 409 for (i=0; i<size; i++) { 410 if (nprocs[2*i+1]) { 411 ierr = MPI_Isend(svalues+starts[i],nprocs[2*i],MPIU_INT,i,tag,comm,send_waits+count++);CHKERRQ(ierr); 412 } 413 } 414 ierr = PetscFree(starts);CHKERRQ(ierr); 415 416 base = owners[rank]; 417 418 /* wait on receives */ 419 ierr = PetscMalloc2(nrecvs,PetscInt,&lens,nrecvs,PetscInt,&source);CHKERRQ(ierr); 420 count = nrecvs; 421 slen = 0; 422 while (count) { 423 ierr = MPI_Waitany(nrecvs,recv_waits,&imdex,&recv_status);CHKERRQ(ierr); 424 /* unpack receives into our local space */ 425 ierr = MPI_Get_count(&recv_status,MPIU_INT,&n);CHKERRQ(ierr); 426 source[imdex] = recv_status.MPI_SOURCE; 427 lens[imdex] = n; 428 slen += n; 429 count--; 430 } 431 ierr = PetscFree(recv_waits);CHKERRQ(ierr); 432 433 /* move the data into the send scatter */ 434 ierr = PetscMalloc((slen+1)*sizeof(PetscInt),&lrows);CHKERRQ(ierr); 435 count = 0; 436 for (i=0; i<nrecvs; i++) { 437 values = rvalues + i*nmax; 438 for (j=0; j<lens[i]; j++) { 439 lrows[count++] = values[j] - base; 440 } 441 } 442 ierr = PetscFree(rvalues);CHKERRQ(ierr); 443 ierr = PetscFree2(lens,source);CHKERRQ(ierr); 444 ierr = PetscFree(owner);CHKERRQ(ierr); 445 ierr = PetscFree(nprocs);CHKERRQ(ierr); 446 447 /* fix right hand side if needed */ 448 if (x && b) { 449 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 450 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 451 for (i=0; i<slen; i++) { 452 bb[lrows[i]] = diag*xx[lrows[i]]; 453 } 454 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 455 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 456 } 457 458 /* actually zap the local rows */ 459 ierr = MatZeroRows(l->A,slen,lrows,0.0,0,0);CHKERRQ(ierr); 460 if (diag != 0.0) { 461 Mat_SeqDense *ll = (Mat_SeqDense*)l->A->data; 462 PetscInt m = ll->lda, i; 463 464 for (i=0; i<slen; i++) { 465 ll->v[lrows[i] + m*(A->cmap->rstart + lrows[i])] = diag; 466 } 467 } 468 ierr = PetscFree(lrows);CHKERRQ(ierr); 469 470 /* wait on sends */ 471 if (nsends) { 472 ierr = PetscMalloc(nsends*sizeof(MPI_Status),&send_status);CHKERRQ(ierr); 473 ierr = MPI_Waitall(nsends,send_waits,send_status);CHKERRQ(ierr); 474 ierr = PetscFree(send_status);CHKERRQ(ierr); 475 } 476 ierr = PetscFree(send_waits);CHKERRQ(ierr); 477 ierr = PetscFree(svalues);CHKERRQ(ierr); 478 479 PetscFunctionReturn(0); 480 } 481 482 #undef __FUNCT__ 483 #define __FUNCT__ "MatMult_MPIDense" 484 PetscErrorCode MatMult_MPIDense(Mat mat,Vec xx,Vec yy) 485 { 486 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 487 PetscErrorCode ierr; 488 489 PetscFunctionBegin; 490 ierr = VecScatterBegin(mdn->Mvctx,xx,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 491 ierr = VecScatterEnd(mdn->Mvctx,xx,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 492 ierr = MatMult_SeqDense(mdn->A,mdn->lvec,yy);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 #undef __FUNCT__ 497 #define __FUNCT__ "MatMultAdd_MPIDense" 498 PetscErrorCode MatMultAdd_MPIDense(Mat mat,Vec xx,Vec yy,Vec zz) 499 { 500 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 501 PetscErrorCode ierr; 502 503 PetscFunctionBegin; 504 ierr = VecScatterBegin(mdn->Mvctx,xx,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 505 ierr = VecScatterEnd(mdn->Mvctx,xx,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 506 ierr = MatMultAdd_SeqDense(mdn->A,mdn->lvec,yy,zz);CHKERRQ(ierr); 507 PetscFunctionReturn(0); 508 } 509 510 #undef __FUNCT__ 511 #define __FUNCT__ "MatMultTranspose_MPIDense" 512 PetscErrorCode MatMultTranspose_MPIDense(Mat A,Vec xx,Vec yy) 513 { 514 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 515 PetscErrorCode ierr; 516 PetscScalar zero = 0.0; 517 518 PetscFunctionBegin; 519 ierr = VecSet(yy,zero);CHKERRQ(ierr); 520 ierr = MatMultTranspose_SeqDense(a->A,xx,a->lvec);CHKERRQ(ierr); 521 ierr = VecScatterBegin(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 522 ierr = VecScatterEnd(a->Mvctx,a->lvec,yy,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 523 PetscFunctionReturn(0); 524 } 525 526 #undef __FUNCT__ 527 #define __FUNCT__ "MatMultTransposeAdd_MPIDense" 528 PetscErrorCode MatMultTransposeAdd_MPIDense(Mat A,Vec xx,Vec yy,Vec zz) 529 { 530 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 531 PetscErrorCode ierr; 532 533 PetscFunctionBegin; 534 ierr = VecCopy(yy,zz);CHKERRQ(ierr); 535 ierr = MatMultTranspose_SeqDense(a->A,xx,a->lvec);CHKERRQ(ierr); 536 ierr = VecScatterBegin(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 537 ierr = VecScatterEnd(a->Mvctx,a->lvec,zz,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 538 PetscFunctionReturn(0); 539 } 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "MatGetDiagonal_MPIDense" 543 PetscErrorCode MatGetDiagonal_MPIDense(Mat A,Vec v) 544 { 545 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 546 Mat_SeqDense *aloc = (Mat_SeqDense*)a->A->data; 547 PetscErrorCode ierr; 548 PetscInt len,i,n,m = A->rmap->n,radd; 549 PetscScalar *x,zero = 0.0; 550 551 PetscFunctionBegin; 552 ierr = VecSet(v,zero);CHKERRQ(ierr); 553 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 554 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 555 if (n != A->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 556 len = PetscMin(a->A->rmap->n,a->A->cmap->n); 557 radd = A->rmap->rstart*m; 558 for (i=0; i<len; i++) { 559 x[i] = aloc->v[radd + i*m + i]; 560 } 561 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "MatDestroy_MPIDense" 567 PetscErrorCode MatDestroy_MPIDense(Mat mat) 568 { 569 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 570 PetscErrorCode ierr; 571 #if defined(PETSC_HAVE_PLAPACK) 572 Mat_Plapack *lu=(Mat_Plapack*)mat->spptr; 573 #endif 574 575 PetscFunctionBegin; 576 577 #if defined(PETSC_USE_LOG) 578 PetscLogObjectState((PetscObject)mat,"Rows=%D, Cols=%D",mat->rmap->N,mat->cmap->N); 579 #endif 580 ierr = MatStashDestroy_Private(&mat->stash);CHKERRQ(ierr); 581 ierr = MatDestroy(&mdn->A);CHKERRQ(ierr); 582 ierr = VecDestroy(&mdn->lvec);CHKERRQ(ierr); 583 ierr = VecScatterDestroy(&mdn->Mvctx);CHKERRQ(ierr); 584 #if defined(PETSC_HAVE_PLAPACK) 585 if (lu) { 586 ierr = PLA_Obj_free(&lu->A);CHKERRQ(ierr); 587 ierr = PLA_Obj_free (&lu->pivots);CHKERRQ(ierr); 588 ierr = PLA_Temp_free(&lu->templ);CHKERRQ(ierr); 589 ierr = ISDestroy(&lu->is_pla);CHKERRQ(ierr); 590 ierr = ISDestroy(&lu->is_petsc);CHKERRQ(ierr); 591 ierr = VecScatterDestroy(&lu->ctx);CHKERRQ(ierr); 592 } 593 ierr = PetscFree(mat->spptr);CHKERRQ(ierr); 594 #endif 595 596 ierr = PetscFree(mat->data);CHKERRQ(ierr); 597 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 598 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetDiagonalBlock_C","",PETSC_NULL);CHKERRQ(ierr); 599 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMPIDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 600 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMult_mpiaij_mpidense_C","",PETSC_NULL);CHKERRQ(ierr); 601 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultSymbolic_mpiaij_mpidense_C","",PETSC_NULL);CHKERRQ(ierr); 602 ierr = PetscObjectComposeFunction((PetscObject)mat,"MatMatMultNumeric_mpiaij_mpidense_C","",PETSC_NULL);CHKERRQ(ierr); 603 PetscFunctionReturn(0); 604 } 605 606 #undef __FUNCT__ 607 #define __FUNCT__ "MatView_MPIDense_Binary" 608 static PetscErrorCode MatView_MPIDense_Binary(Mat mat,PetscViewer viewer) 609 { 610 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 611 PetscErrorCode ierr; 612 PetscViewerFormat format; 613 int fd; 614 PetscInt header[4],mmax,N = mat->cmap->N,i,j,m,k; 615 PetscMPIInt rank,tag = ((PetscObject)viewer)->tag,size; 616 PetscScalar *work,*v,*vv; 617 Mat_SeqDense *a = (Mat_SeqDense*)mdn->A->data; 618 619 PetscFunctionBegin; 620 if (mdn->size == 1) { 621 ierr = MatView(mdn->A,viewer);CHKERRQ(ierr); 622 } else { 623 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 624 ierr = MPI_Comm_rank(((PetscObject)mat)->comm,&rank);CHKERRQ(ierr); 625 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&size);CHKERRQ(ierr); 626 627 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 628 if (format == PETSC_VIEWER_NATIVE) { 629 630 if (!rank) { 631 /* store the matrix as a dense matrix */ 632 header[0] = MAT_FILE_CLASSID; 633 header[1] = mat->rmap->N; 634 header[2] = N; 635 header[3] = MATRIX_BINARY_FORMAT_DENSE; 636 ierr = PetscBinaryWrite(fd,header,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 637 638 /* get largest work array needed for transposing array */ 639 mmax = mat->rmap->n; 640 for (i=1; i<size; i++) { 641 mmax = PetscMax(mmax,mat->rmap->range[i+1] - mat->rmap->range[i]); 642 } 643 ierr = PetscMalloc(mmax*N*sizeof(PetscScalar),&work);CHKERRQ(ierr); 644 645 /* write out local array, by rows */ 646 m = mat->rmap->n; 647 v = a->v; 648 for (j=0; j<N; j++) { 649 for (i=0; i<m; i++) { 650 work[j + i*N] = *v++; 651 } 652 } 653 ierr = PetscBinaryWrite(fd,work,m*N,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 654 /* get largest work array to receive messages from other processes, excludes process zero */ 655 mmax = 0; 656 for (i=1; i<size; i++) { 657 mmax = PetscMax(mmax,mat->rmap->range[i+1] - mat->rmap->range[i]); 658 } 659 ierr = PetscMalloc(mmax*N*sizeof(PetscScalar),&vv);CHKERRQ(ierr); 660 for(k = 1; k < size; k++) { 661 v = vv; 662 m = mat->rmap->range[k+1] - mat->rmap->range[k]; 663 ierr = MPILong_Recv(v,m*N,MPIU_SCALAR,k,tag,((PetscObject)mat)->comm);CHKERRQ(ierr); 664 665 for(j = 0; j < N; j++) { 666 for(i = 0; i < m; i++) { 667 work[j + i*N] = *v++; 668 } 669 } 670 ierr = PetscBinaryWrite(fd,work,m*N,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 671 } 672 ierr = PetscFree(work);CHKERRQ(ierr); 673 ierr = PetscFree(vv);CHKERRQ(ierr); 674 } else { 675 ierr = MPILong_Send(a->v,mat->rmap->n*mat->cmap->N,MPIU_SCALAR,0,tag,((PetscObject)mat)->comm);CHKERRQ(ierr); 676 } 677 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"To store a parallel dense matrix you must first call PetscViewerSetFormat(viewer,PETSC_VIEWER_NATIVE)"); 678 } 679 PetscFunctionReturn(0); 680 } 681 682 #undef __FUNCT__ 683 #define __FUNCT__ "MatView_MPIDense_ASCIIorDraworSocket" 684 static PetscErrorCode MatView_MPIDense_ASCIIorDraworSocket(Mat mat,PetscViewer viewer) 685 { 686 Mat_MPIDense *mdn = (Mat_MPIDense*)mat->data; 687 PetscErrorCode ierr; 688 PetscMPIInt size = mdn->size,rank = mdn->rank; 689 const PetscViewerType vtype; 690 PetscBool iascii,isdraw; 691 PetscViewer sviewer; 692 PetscViewerFormat format; 693 #if defined(PETSC_HAVE_PLAPACK) 694 Mat_Plapack *lu; 695 #endif 696 697 PetscFunctionBegin; 698 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 699 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 700 if (iascii) { 701 ierr = PetscViewerGetType(viewer,&vtype);CHKERRQ(ierr); 702 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 703 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 704 MatInfo info; 705 ierr = MatGetInfo(mat,MAT_LOCAL,&info);CHKERRQ(ierr); 706 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);CHKERRQ(ierr); 707 ierr = PetscViewerASCIISynchronizedPrintf(viewer," [%d] local rows %D nz %D nz alloced %D mem %D \n",rank,mat->rmap->n,(PetscInt)info.nz_used,(PetscInt)info.nz_allocated,(PetscInt)info.memory);CHKERRQ(ierr); 708 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 709 ierr = PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);CHKERRQ(ierr); 710 #if defined(PETSC_HAVE_PLAPACK) 711 ierr = PetscViewerASCIIPrintf(viewer,"PLAPACK run parameters:\n");CHKERRQ(ierr); 712 ierr = PetscViewerASCIIPrintf(viewer," Processor mesh: nprows %d, npcols %d\n",Plapack_nprows, Plapack_npcols);CHKERRQ(ierr); 713 ierr = PetscViewerASCIIPrintf(viewer," Error checking: %d\n",Plapack_ierror);CHKERRQ(ierr); 714 ierr = PetscViewerASCIIPrintf(viewer," Algorithmic block size: %d\n",Plapack_nb_alg);CHKERRQ(ierr); 715 if (mat->factortype){ 716 lu=(Mat_Plapack*)(mat->spptr); 717 ierr = PetscViewerASCIIPrintf(viewer," Distr. block size nb: %d \n",lu->nb);CHKERRQ(ierr); 718 } 719 #else 720 ierr = VecScatterView(mdn->Mvctx,viewer);CHKERRQ(ierr); 721 #endif 722 PetscFunctionReturn(0); 723 } else if (format == PETSC_VIEWER_ASCII_INFO) { 724 PetscFunctionReturn(0); 725 } 726 } else if (isdraw) { 727 PetscDraw draw; 728 PetscBool isnull; 729 730 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 731 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 732 if (isnull) PetscFunctionReturn(0); 733 } 734 735 if (size == 1) { 736 ierr = MatView(mdn->A,viewer);CHKERRQ(ierr); 737 } else { 738 /* assemble the entire matrix onto first processor. */ 739 Mat A; 740 PetscInt M = mat->rmap->N,N = mat->cmap->N,m,row,i,nz; 741 PetscInt *cols; 742 PetscScalar *vals; 743 744 ierr = MatCreate(((PetscObject)mat)->comm,&A);CHKERRQ(ierr); 745 if (!rank) { 746 ierr = MatSetSizes(A,M,N,M,N);CHKERRQ(ierr); 747 } else { 748 ierr = MatSetSizes(A,0,0,M,N);CHKERRQ(ierr); 749 } 750 /* Since this is a temporary matrix, MATMPIDENSE instead of ((PetscObject)A)->type_name here is probably acceptable. */ 751 ierr = MatSetType(A,MATMPIDENSE);CHKERRQ(ierr); 752 ierr = MatMPIDenseSetPreallocation(A,PETSC_NULL); 753 ierr = PetscLogObjectParent(mat,A);CHKERRQ(ierr); 754 755 /* Copy the matrix ... This isn't the most efficient means, 756 but it's quick for now */ 757 A->insertmode = INSERT_VALUES; 758 row = mat->rmap->rstart; m = mdn->A->rmap->n; 759 for (i=0; i<m; i++) { 760 ierr = MatGetRow_MPIDense(mat,row,&nz,&cols,&vals);CHKERRQ(ierr); 761 ierr = MatSetValues_MPIDense(A,1,&row,nz,cols,vals,INSERT_VALUES);CHKERRQ(ierr); 762 ierr = MatRestoreRow_MPIDense(mat,row,&nz,&cols,&vals);CHKERRQ(ierr); 763 row++; 764 } 765 766 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 767 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 768 ierr = PetscViewerGetSingleton(viewer,&sviewer);CHKERRQ(ierr); 769 if (!rank) { 770 ierr = PetscObjectSetName((PetscObject)((Mat_MPIDense*)(A->data))->A,((PetscObject)mat)->name);CHKERRQ(ierr); 771 /* Set the type name to MATMPIDense so that the correct type can be printed out by PetscObjectPrintClassNamePrefixType() in MatView_SeqDense_ASCII()*/ 772 PetscStrcpy(((PetscObject)((Mat_MPIDense*)(A->data))->A)->type_name,MATMPIDENSE); 773 ierr = MatView(((Mat_MPIDense*)(A->data))->A,sviewer);CHKERRQ(ierr); 774 } 775 ierr = PetscViewerRestoreSingleton(viewer,&sviewer);CHKERRQ(ierr); 776 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 777 ierr = MatDestroy(&A);CHKERRQ(ierr); 778 } 779 PetscFunctionReturn(0); 780 } 781 782 #undef __FUNCT__ 783 #define __FUNCT__ "MatView_MPIDense" 784 PetscErrorCode MatView_MPIDense(Mat mat,PetscViewer viewer) 785 { 786 PetscErrorCode ierr; 787 PetscBool iascii,isbinary,isdraw,issocket; 788 789 PetscFunctionBegin; 790 791 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 792 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 793 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);CHKERRQ(ierr); 794 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 795 796 if (iascii || issocket || isdraw) { 797 ierr = MatView_MPIDense_ASCIIorDraworSocket(mat,viewer);CHKERRQ(ierr); 798 } else if (isbinary) { 799 ierr = MatView_MPIDense_Binary(mat,viewer);CHKERRQ(ierr); 800 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by MPI dense matrix",((PetscObject)viewer)->type_name); 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNCT__ 805 #define __FUNCT__ "MatGetInfo_MPIDense" 806 PetscErrorCode MatGetInfo_MPIDense(Mat A,MatInfoType flag,MatInfo *info) 807 { 808 Mat_MPIDense *mat = (Mat_MPIDense*)A->data; 809 Mat mdn = mat->A; 810 PetscErrorCode ierr; 811 PetscReal isend[5],irecv[5]; 812 813 PetscFunctionBegin; 814 info->block_size = 1.0; 815 ierr = MatGetInfo(mdn,MAT_LOCAL,info);CHKERRQ(ierr); 816 isend[0] = info->nz_used; isend[1] = info->nz_allocated; isend[2] = info->nz_unneeded; 817 isend[3] = info->memory; isend[4] = info->mallocs; 818 if (flag == MAT_LOCAL) { 819 info->nz_used = isend[0]; 820 info->nz_allocated = isend[1]; 821 info->nz_unneeded = isend[2]; 822 info->memory = isend[3]; 823 info->mallocs = isend[4]; 824 } else if (flag == MAT_GLOBAL_MAX) { 825 ierr = MPI_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_MAX,((PetscObject)A)->comm);CHKERRQ(ierr); 826 info->nz_used = irecv[0]; 827 info->nz_allocated = irecv[1]; 828 info->nz_unneeded = irecv[2]; 829 info->memory = irecv[3]; 830 info->mallocs = irecv[4]; 831 } else if (flag == MAT_GLOBAL_SUM) { 832 ierr = MPI_Allreduce(isend,irecv,5,MPIU_REAL,MPIU_SUM,((PetscObject)A)->comm);CHKERRQ(ierr); 833 info->nz_used = irecv[0]; 834 info->nz_allocated = irecv[1]; 835 info->nz_unneeded = irecv[2]; 836 info->memory = irecv[3]; 837 info->mallocs = irecv[4]; 838 } 839 info->fill_ratio_given = 0; /* no parallel LU/ILU/Cholesky */ 840 info->fill_ratio_needed = 0; 841 info->factor_mallocs = 0; 842 PetscFunctionReturn(0); 843 } 844 845 #undef __FUNCT__ 846 #define __FUNCT__ "MatSetOption_MPIDense" 847 PetscErrorCode MatSetOption_MPIDense(Mat A,MatOption op,PetscBool flg) 848 { 849 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 850 PetscErrorCode ierr; 851 852 PetscFunctionBegin; 853 switch (op) { 854 case MAT_NEW_NONZERO_LOCATIONS: 855 case MAT_NEW_NONZERO_LOCATION_ERR: 856 case MAT_NEW_NONZERO_ALLOCATION_ERR: 857 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 858 break; 859 case MAT_ROW_ORIENTED: 860 a->roworiented = flg; 861 ierr = MatSetOption(a->A,op,flg);CHKERRQ(ierr); 862 break; 863 case MAT_NEW_DIAGONALS: 864 case MAT_USE_HASH_TABLE: 865 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 866 break; 867 case MAT_IGNORE_OFF_PROC_ENTRIES: 868 a->donotstash = flg; 869 break; 870 case MAT_SYMMETRIC: 871 case MAT_STRUCTURALLY_SYMMETRIC: 872 case MAT_HERMITIAN: 873 case MAT_SYMMETRY_ETERNAL: 874 case MAT_IGNORE_LOWER_TRIANGULAR: 875 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 876 break; 877 default: 878 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 879 } 880 PetscFunctionReturn(0); 881 } 882 883 884 #undef __FUNCT__ 885 #define __FUNCT__ "MatDiagonalScale_MPIDense" 886 PetscErrorCode MatDiagonalScale_MPIDense(Mat A,Vec ll,Vec rr) 887 { 888 Mat_MPIDense *mdn = (Mat_MPIDense*)A->data; 889 Mat_SeqDense *mat = (Mat_SeqDense*)mdn->A->data; 890 PetscScalar *l,*r,x,*v; 891 PetscErrorCode ierr; 892 PetscInt i,j,s2a,s3a,s2,s3,m=mdn->A->rmap->n,n=mdn->A->cmap->n; 893 894 PetscFunctionBegin; 895 ierr = MatGetLocalSize(A,&s2,&s3);CHKERRQ(ierr); 896 if (ll) { 897 ierr = VecGetLocalSize(ll,&s2a);CHKERRQ(ierr); 898 if (s2a != s2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector non-conforming local size, %d != %d.", s2a, s2); 899 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 900 for (i=0; i<m; i++) { 901 x = l[i]; 902 v = mat->v + i; 903 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 904 } 905 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 906 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 907 } 908 if (rr) { 909 ierr = VecGetLocalSize(rr,&s3a);CHKERRQ(ierr); 910 if (s3a != s3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec non-conforming local size, %d != %d.", s3a, s3); 911 ierr = VecScatterBegin(mdn->Mvctx,rr,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 912 ierr = VecScatterEnd(mdn->Mvctx,rr,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 913 ierr = VecGetArray(mdn->lvec,&r);CHKERRQ(ierr); 914 for (i=0; i<n; i++) { 915 x = r[i]; 916 v = mat->v + i*m; 917 for (j=0; j<m; j++) { (*v++) *= x;} 918 } 919 ierr = VecRestoreArray(mdn->lvec,&r);CHKERRQ(ierr); 920 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 921 } 922 PetscFunctionReturn(0); 923 } 924 925 #undef __FUNCT__ 926 #define __FUNCT__ "MatNorm_MPIDense" 927 PetscErrorCode MatNorm_MPIDense(Mat A,NormType type,PetscReal *nrm) 928 { 929 Mat_MPIDense *mdn = (Mat_MPIDense*)A->data; 930 Mat_SeqDense *mat = (Mat_SeqDense*)mdn->A->data; 931 PetscErrorCode ierr; 932 PetscInt i,j; 933 PetscReal sum = 0.0; 934 PetscScalar *v = mat->v; 935 936 PetscFunctionBegin; 937 if (mdn->size == 1) { 938 ierr = MatNorm(mdn->A,type,nrm);CHKERRQ(ierr); 939 } else { 940 if (type == NORM_FROBENIUS) { 941 for (i=0; i<mdn->A->cmap->n*mdn->A->rmap->n; i++) { 942 #if defined(PETSC_USE_COMPLEX) 943 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 944 #else 945 sum += (*v)*(*v); v++; 946 #endif 947 } 948 ierr = MPI_Allreduce(&sum,nrm,1,MPIU_REAL,MPIU_SUM,((PetscObject)A)->comm);CHKERRQ(ierr); 949 *nrm = PetscSqrtReal(*nrm); 950 ierr = PetscLogFlops(2.0*mdn->A->cmap->n*mdn->A->rmap->n);CHKERRQ(ierr); 951 } else if (type == NORM_1) { 952 PetscReal *tmp,*tmp2; 953 ierr = PetscMalloc2(A->cmap->N,PetscReal,&tmp,A->cmap->N,PetscReal,&tmp2);CHKERRQ(ierr); 954 ierr = PetscMemzero(tmp,A->cmap->N*sizeof(PetscReal));CHKERRQ(ierr); 955 ierr = PetscMemzero(tmp2,A->cmap->N*sizeof(PetscReal));CHKERRQ(ierr); 956 *nrm = 0.0; 957 v = mat->v; 958 for (j=0; j<mdn->A->cmap->n; j++) { 959 for (i=0; i<mdn->A->rmap->n; i++) { 960 tmp[j] += PetscAbsScalar(*v); v++; 961 } 962 } 963 ierr = MPI_Allreduce(tmp,tmp2,A->cmap->N,MPIU_REAL,MPIU_SUM,((PetscObject)A)->comm);CHKERRQ(ierr); 964 for (j=0; j<A->cmap->N; j++) { 965 if (tmp2[j] > *nrm) *nrm = tmp2[j]; 966 } 967 ierr = PetscFree2(tmp,tmp);CHKERRQ(ierr); 968 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 969 } else if (type == NORM_INFINITY) { /* max row norm */ 970 PetscReal ntemp; 971 ierr = MatNorm(mdn->A,type,&ntemp);CHKERRQ(ierr); 972 ierr = MPI_Allreduce(&ntemp,nrm,1,MPIU_REAL,MPIU_MAX,((PetscObject)A)->comm);CHKERRQ(ierr); 973 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support for two norm"); 974 } 975 PetscFunctionReturn(0); 976 } 977 978 #undef __FUNCT__ 979 #define __FUNCT__ "MatTranspose_MPIDense" 980 PetscErrorCode MatTranspose_MPIDense(Mat A,MatReuse reuse,Mat *matout) 981 { 982 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 983 Mat_SeqDense *Aloc = (Mat_SeqDense*)a->A->data; 984 Mat B; 985 PetscInt M = A->rmap->N,N = A->cmap->N,m,n,*rwork,rstart = A->rmap->rstart; 986 PetscErrorCode ierr; 987 PetscInt j,i; 988 PetscScalar *v; 989 990 PetscFunctionBegin; 991 if (reuse == MAT_REUSE_MATRIX && A == *matout && M != N) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Supports square matrix only in-place"); 992 if (reuse == MAT_INITIAL_MATRIX || A == *matout) { 993 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 994 ierr = MatSetSizes(B,A->cmap->n,A->rmap->n,N,M);CHKERRQ(ierr); 995 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 996 ierr = MatMPIDenseSetPreallocation(B,PETSC_NULL);CHKERRQ(ierr); 997 } else { 998 B = *matout; 999 } 1000 1001 m = a->A->rmap->n; n = a->A->cmap->n; v = Aloc->v; 1002 ierr = PetscMalloc(m*sizeof(PetscInt),&rwork);CHKERRQ(ierr); 1003 for (i=0; i<m; i++) rwork[i] = rstart + i; 1004 for (j=0; j<n; j++) { 1005 ierr = MatSetValues(B,1,&j,m,rwork,v,INSERT_VALUES);CHKERRQ(ierr); 1006 v += m; 1007 } 1008 ierr = PetscFree(rwork);CHKERRQ(ierr); 1009 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1010 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 if (reuse == MAT_INITIAL_MATRIX || *matout != A) { 1012 *matout = B; 1013 } else { 1014 ierr = MatHeaderMerge(A,B);CHKERRQ(ierr); 1015 } 1016 PetscFunctionReturn(0); 1017 } 1018 1019 1020 static PetscErrorCode MatDuplicate_MPIDense(Mat,MatDuplicateOption,Mat *); 1021 extern PetscErrorCode MatScale_MPIDense(Mat,PetscScalar); 1022 1023 #undef __FUNCT__ 1024 #define __FUNCT__ "MatSetUpPreallocation_MPIDense" 1025 PetscErrorCode MatSetUpPreallocation_MPIDense(Mat A) 1026 { 1027 PetscErrorCode ierr; 1028 1029 PetscFunctionBegin; 1030 ierr = MatMPIDenseSetPreallocation(A,0);CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 1034 #if defined(PETSC_HAVE_PLAPACK) 1035 1036 #undef __FUNCT__ 1037 #define __FUNCT__ "MatMPIDenseCopyToPlapack" 1038 PetscErrorCode MatMPIDenseCopyToPlapack(Mat A,Mat F) 1039 { 1040 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1041 PetscErrorCode ierr; 1042 PetscInt M=A->cmap->N,m=A->rmap->n,rstart; 1043 PetscScalar *array; 1044 PetscReal one = 1.0; 1045 1046 PetscFunctionBegin; 1047 /* Copy A into F->lu->A */ 1048 ierr = PLA_Obj_set_to_zero(lu->A);CHKERRQ(ierr); 1049 ierr = PLA_API_begin();CHKERRQ(ierr); 1050 ierr = PLA_Obj_API_open(lu->A);CHKERRQ(ierr); 1051 ierr = MatGetOwnershipRange(A,&rstart,PETSC_NULL);CHKERRQ(ierr); 1052 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1053 ierr = PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0);CHKERRQ(ierr); 1054 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1055 ierr = PLA_Obj_API_close(lu->A);CHKERRQ(ierr); 1056 ierr = PLA_API_end();CHKERRQ(ierr); 1057 lu->rstart = rstart; 1058 PetscFunctionReturn(0); 1059 } 1060 1061 #undef __FUNCT__ 1062 #define __FUNCT__ "MatMPIDenseCopyFromPlapack" 1063 PetscErrorCode MatMPIDenseCopyFromPlapack(Mat F,Mat A) 1064 { 1065 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1066 PetscErrorCode ierr; 1067 PetscInt M=A->cmap->N,m=A->rmap->n,rstart; 1068 PetscScalar *array; 1069 PetscReal one = 1.0; 1070 1071 PetscFunctionBegin; 1072 /* Copy F into A->lu->A */ 1073 ierr = MatZeroEntries(A);CHKERRQ(ierr); 1074 ierr = PLA_API_begin();CHKERRQ(ierr); 1075 ierr = PLA_Obj_API_open(lu->A);CHKERRQ(ierr); 1076 ierr = MatGetOwnershipRange(A,&rstart,PETSC_NULL);CHKERRQ(ierr); 1077 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1078 ierr = PLA_API_axpy_global_to_matrix(m,M, &one,lu->A,rstart,0,(void *)array,m);CHKERRQ(ierr); 1079 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1080 ierr = PLA_Obj_API_close(lu->A);CHKERRQ(ierr); 1081 ierr = PLA_API_end();CHKERRQ(ierr); 1082 lu->rstart = rstart; 1083 PetscFunctionReturn(0); 1084 } 1085 1086 #undef __FUNCT__ 1087 #define __FUNCT__ "MatMatMultNumeric_MPIDense_MPIDense" 1088 PetscErrorCode MatMatMultNumeric_MPIDense_MPIDense(Mat A,Mat B,Mat C) 1089 { 1090 PetscErrorCode ierr; 1091 Mat_Plapack *luA = (Mat_Plapack*)A->spptr; 1092 Mat_Plapack *luB = (Mat_Plapack*)B->spptr; 1093 Mat_Plapack *luC = (Mat_Plapack*)C->spptr; 1094 PLA_Obj alpha = NULL,beta = NULL; 1095 1096 PetscFunctionBegin; 1097 ierr = MatMPIDenseCopyToPlapack(A,A);CHKERRQ(ierr); 1098 ierr = MatMPIDenseCopyToPlapack(B,B);CHKERRQ(ierr); 1099 1100 /* 1101 ierr = PLA_Global_show("A = ",luA->A,"%g ","");CHKERRQ(ierr); 1102 ierr = PLA_Global_show("B = ",luB->A,"%g ","");CHKERRQ(ierr); 1103 */ 1104 1105 /* do the multiply in PLA */ 1106 ierr = PLA_Create_constants_conf_to(luA->A,NULL,NULL,&alpha);CHKERRQ(ierr); 1107 ierr = PLA_Create_constants_conf_to(luC->A,NULL,&beta,NULL);CHKERRQ(ierr); 1108 CHKMEMQ; 1109 1110 ierr = PLA_Gemm(PLA_NO_TRANSPOSE,PLA_NO_TRANSPOSE,alpha,luA->A,luB->A,beta,luC->A); /* CHKERRQ(ierr); */ 1111 CHKMEMQ; 1112 ierr = PLA_Obj_free(&alpha);CHKERRQ(ierr); 1113 ierr = PLA_Obj_free(&beta);CHKERRQ(ierr); 1114 1115 /* 1116 ierr = PLA_Global_show("C = ",luC->A,"%g ","");CHKERRQ(ierr); 1117 */ 1118 ierr = MatMPIDenseCopyFromPlapack(C,C);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "MatMatMultSymbolic_MPIDense_MPIDense" 1124 PetscErrorCode MatMatMultSymbolic_MPIDense_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C) 1125 { 1126 PetscErrorCode ierr; 1127 PetscInt m=A->rmap->n,n=B->cmap->n; 1128 Mat Cmat; 1129 1130 PetscFunctionBegin; 1131 if (A->cmap->n != B->rmap->n) SETERRQ2(((PetscObject)A)->comm,PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 1132 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_LIB,"Due to apparent bugs in PLAPACK,this is not currently supported"); 1133 ierr = MatCreate(((PetscObject)B)->comm,&Cmat);CHKERRQ(ierr); 1134 ierr = MatSetSizes(Cmat,m,n,A->rmap->N,B->cmap->N);CHKERRQ(ierr); 1135 ierr = MatSetType(Cmat,MATMPIDENSE);CHKERRQ(ierr); 1136 ierr = MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1137 ierr = MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1138 1139 *C = Cmat; 1140 PetscFunctionReturn(0); 1141 } 1142 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "MatMatMult_MPIDense_MPIDense" 1145 PetscErrorCode MatMatMult_MPIDense_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1146 { 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 if (scall == MAT_INITIAL_MATRIX){ 1151 ierr = MatMatMultSymbolic_MPIDense_MPIDense(A,B,fill,C);CHKERRQ(ierr); 1152 } 1153 ierr = MatMatMultNumeric_MPIDense_MPIDense(A,B,*C);CHKERRQ(ierr); 1154 PetscFunctionReturn(0); 1155 } 1156 1157 #undef __FUNCT__ 1158 #define __FUNCT__ "MatSolve_MPIDense" 1159 PetscErrorCode MatSolve_MPIDense(Mat A,Vec b,Vec x) 1160 { 1161 MPI_Comm comm = ((PetscObject)A)->comm; 1162 Mat_Plapack *lu = (Mat_Plapack*)A->spptr; 1163 PetscErrorCode ierr; 1164 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,i,j,*idx_pla,*idx_petsc,loc_m,loc_stride; 1165 PetscScalar *array; 1166 PetscReal one = 1.0; 1167 PetscMPIInt size,rank,r_rank,r_nproc,c_rank,c_nproc;; 1168 PLA_Obj v_pla = NULL; 1169 PetscScalar *loc_buf; 1170 Vec loc_x; 1171 1172 PetscFunctionBegin; 1173 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1174 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1175 1176 /* Create PLAPACK vector objects, then copy b into PLAPACK b */ 1177 PLA_Mvector_create(lu->datatype,M,1,lu->templ,PLA_ALIGN_FIRST,&v_pla); 1178 PLA_Obj_set_to_zero(v_pla); 1179 1180 /* Copy b into rhs_pla */ 1181 PLA_API_begin(); 1182 PLA_Obj_API_open(v_pla); 1183 ierr = VecGetArray(b,&array);CHKERRQ(ierr); 1184 PLA_API_axpy_vector_to_global(m,&one,(void *)array,1,v_pla,lu->rstart); 1185 ierr = VecRestoreArray(b,&array);CHKERRQ(ierr); 1186 PLA_Obj_API_close(v_pla); 1187 PLA_API_end(); 1188 1189 if (A->factortype == MAT_FACTOR_LU){ 1190 /* Apply the permutations to the right hand sides */ 1191 PLA_Apply_pivots_to_rows (v_pla,lu->pivots); 1192 1193 /* Solve L y = b, overwriting b with y */ 1194 PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_UNIT_DIAG,lu->A,v_pla ); 1195 1196 /* Solve U x = y (=b), overwriting b with x */ 1197 PLA_Trsv( PLA_UPPER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla ); 1198 } else { /* MAT_FACTOR_CHOLESKY */ 1199 PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla); 1200 PLA_Trsv( PLA_LOWER_TRIANGULAR,(lu->datatype == MPI_DOUBLE ? PLA_TRANSPOSE : PLA_CONJUGATE_TRANSPOSE), 1201 PLA_NONUNIT_DIAG,lu->A,v_pla); 1202 } 1203 1204 /* Copy PLAPACK x into Petsc vector x */ 1205 PLA_Obj_local_length(v_pla, &loc_m); 1206 PLA_Obj_local_buffer(v_pla, (void**)&loc_buf); 1207 PLA_Obj_local_stride(v_pla, &loc_stride); 1208 /* 1209 PetscPrintf(PETSC_COMM_SELF," [%d] b - local_m %d local_stride %d, loc_buf: %g %g, nb: %d\n",rank,loc_m,loc_stride,loc_buf[0],loc_buf[(loc_m-1)*loc_stride],lu->nb); 1210 */ 1211 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,loc_m*loc_stride,loc_buf,&loc_x);CHKERRQ(ierr); 1212 if (!lu->pla_solved){ 1213 1214 PLA_Temp_comm_row_info(lu->templ,&Plapack_comm_2d,&r_rank,&r_nproc); 1215 PLA_Temp_comm_col_info(lu->templ,&Plapack_comm_2d,&c_rank,&c_nproc); 1216 1217 /* Create IS and cts for VecScatterring */ 1218 PLA_Obj_local_length(v_pla, &loc_m); 1219 PLA_Obj_local_stride(v_pla, &loc_stride); 1220 ierr = PetscMalloc2(loc_m,PetscInt,&idx_pla,loc_m,PetscInt,&idx_petsc);CHKERRQ(ierr); 1221 1222 rstart = (r_rank*c_nproc+c_rank)*lu->nb; 1223 for (i=0; i<loc_m; i+=lu->nb){ 1224 j = 0; 1225 while (j < lu->nb && i+j < loc_m){ 1226 idx_petsc[i+j] = rstart + j; j++; 1227 } 1228 rstart += size*lu->nb; 1229 } 1230 1231 for (i=0; i<loc_m; i++) idx_pla[i] = i*loc_stride; 1232 1233 ierr = ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_pla,PETSC_COPY_VALUES,&lu->is_pla);CHKERRQ(ierr); 1234 ierr = ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_petsc,PETSC_COPY_VALUES,&lu->is_petsc);CHKERRQ(ierr); 1235 ierr = PetscFree2(idx_pla,idx_petsc);CHKERRQ(ierr); 1236 ierr = VecScatterCreate(loc_x,lu->is_pla,x,lu->is_petsc,&lu->ctx);CHKERRQ(ierr); 1237 } 1238 ierr = VecScatterBegin(lu->ctx,loc_x,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1239 ierr = VecScatterEnd(lu->ctx,loc_x,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1240 1241 /* Free data */ 1242 ierr = VecDestroy(&loc_x);CHKERRQ(ierr); 1243 PLA_Obj_free(&v_pla); 1244 1245 lu->pla_solved = PETSC_TRUE; 1246 PetscFunctionReturn(0); 1247 } 1248 1249 #undef __FUNCT__ 1250 #define __FUNCT__ "MatLUFactorNumeric_MPIDense" 1251 PetscErrorCode MatLUFactorNumeric_MPIDense(Mat F,Mat A,const MatFactorInfo *info) 1252 { 1253 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1254 PetscErrorCode ierr; 1255 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,rend; 1256 PetscInt info_pla=0; 1257 PetscScalar *array,one = 1.0; 1258 1259 PetscFunctionBegin; 1260 if (lu->mstruct == SAME_NONZERO_PATTERN){ 1261 PLA_Obj_free(&lu->A); 1262 PLA_Obj_free (&lu->pivots); 1263 } 1264 /* Create PLAPACK matrix object */ 1265 lu->A = NULL; lu->pivots = NULL; 1266 PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A); 1267 PLA_Obj_set_to_zero(lu->A); 1268 PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots); 1269 1270 /* Copy A into lu->A */ 1271 PLA_API_begin(); 1272 PLA_Obj_API_open(lu->A); 1273 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1274 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1275 PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0); 1276 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1277 PLA_Obj_API_close(lu->A); 1278 PLA_API_end(); 1279 1280 /* Factor P A -> L U overwriting lower triangular portion of A with L, upper, U */ 1281 info_pla = PLA_LU(lu->A,lu->pivots); 1282 if (info_pla != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot encountered at row %d from PLA_LU()",info_pla); 1283 1284 lu->rstart = rstart; 1285 lu->mstruct = SAME_NONZERO_PATTERN; 1286 F->ops->solve = MatSolve_MPIDense; 1287 F->assembled = PETSC_TRUE; /* required by -ksp_view */ 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "MatCholeskyFactorNumeric_MPIDense" 1293 PetscErrorCode MatCholeskyFactorNumeric_MPIDense(Mat F,Mat A,const MatFactorInfo *info) 1294 { 1295 Mat_Plapack *lu = (Mat_Plapack*)F->spptr; 1296 PetscErrorCode ierr; 1297 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,rend; 1298 PetscInt info_pla=0; 1299 PetscScalar *array,one = 1.0; 1300 1301 PetscFunctionBegin; 1302 if (lu->mstruct == SAME_NONZERO_PATTERN){ 1303 PLA_Obj_free(&lu->A); 1304 } 1305 /* Create PLAPACK matrix object */ 1306 lu->A = NULL; 1307 lu->pivots = NULL; 1308 PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A); 1309 1310 /* Copy A into lu->A */ 1311 PLA_API_begin(); 1312 PLA_Obj_API_open(lu->A); 1313 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1314 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1315 PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0); 1316 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1317 PLA_Obj_API_close(lu->A); 1318 PLA_API_end(); 1319 1320 /* Factor P A -> Chol */ 1321 info_pla = PLA_Chol(PLA_LOWER_TRIANGULAR,lu->A); 1322 if (info_pla != 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_MAT_CH_ZRPVT,"Nonpositive definite matrix detected at row %d from PLA_Chol()",info_pla); 1323 1324 lu->rstart = rstart; 1325 lu->mstruct = SAME_NONZERO_PATTERN; 1326 F->ops->solve = MatSolve_MPIDense; 1327 F->assembled = PETSC_TRUE; /* required by -ksp_view */ 1328 PetscFunctionReturn(0); 1329 } 1330 1331 /* Note the Petsc perm permutation is ignored */ 1332 #undef __FUNCT__ 1333 #define __FUNCT__ "MatCholeskyFactorSymbolic_MPIDense" 1334 PetscErrorCode MatCholeskyFactorSymbolic_MPIDense(Mat F,Mat A,IS perm,const MatFactorInfo *info) 1335 { 1336 PetscErrorCode ierr; 1337 PetscBool issymmetric,set; 1338 1339 PetscFunctionBegin; 1340 ierr = MatIsSymmetricKnown(A,&set,&issymmetric);CHKERRQ(ierr); 1341 if (!set || !issymmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be set as MAT_SYMMETRIC for CholeskyFactor()"); 1342 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_MPIDense; 1343 PetscFunctionReturn(0); 1344 } 1345 1346 /* Note the Petsc r and c permutations are ignored */ 1347 #undef __FUNCT__ 1348 #define __FUNCT__ "MatLUFactorSymbolic_MPIDense" 1349 PetscErrorCode MatLUFactorSymbolic_MPIDense(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1350 { 1351 PetscErrorCode ierr; 1352 PetscInt M = A->rmap->N; 1353 Mat_Plapack *lu; 1354 1355 PetscFunctionBegin; 1356 lu = (Mat_Plapack*)F->spptr; 1357 ierr = PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots);CHKERRQ(ierr); 1358 F->ops->lufactornumeric = MatLUFactorNumeric_MPIDense; 1359 PetscFunctionReturn(0); 1360 } 1361 1362 EXTERN_C_BEGIN 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "MatFactorGetSolverPackage_mpidense_plapack" 1365 PetscErrorCode MatFactorGetSolverPackage_mpidense_plapack(Mat A,const MatSolverPackage *type) 1366 { 1367 PetscFunctionBegin; 1368 *type = MATSOLVERPLAPACK; 1369 PetscFunctionReturn(0); 1370 } 1371 EXTERN_C_END 1372 1373 EXTERN_C_BEGIN 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "MatGetFactor_mpidense_plapack" 1376 PetscErrorCode MatGetFactor_mpidense_plapack(Mat A,MatFactorType ftype,Mat *F) 1377 { 1378 PetscErrorCode ierr; 1379 Mat_Plapack *lu; 1380 PetscMPIInt size; 1381 PetscInt M=A->rmap->N; 1382 1383 PetscFunctionBegin; 1384 /* Create the factorization matrix */ 1385 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 1386 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1387 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 1388 ierr = PetscNewLog(*F,Mat_Plapack,&lu);CHKERRQ(ierr); 1389 (*F)->spptr = (void*)lu; 1390 1391 /* Set default Plapack parameters */ 1392 ierr = MPI_Comm_size(((PetscObject)A)->comm,&size);CHKERRQ(ierr); 1393 lu->nb = M/size; 1394 if (M - lu->nb*size) lu->nb++; /* without cyclic distribution */ 1395 1396 /* Set runtime options */ 1397 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"PLAPACK Options","Mat");CHKERRQ(ierr); 1398 ierr = PetscOptionsInt("-mat_plapack_nb","block size of template vector","None",lu->nb,&lu->nb,PETSC_NULL);CHKERRQ(ierr); 1399 PetscOptionsEnd(); 1400 1401 /* Create object distribution template */ 1402 lu->templ = NULL; 1403 ierr = PLA_Temp_create(lu->nb, 0, &lu->templ);CHKERRQ(ierr); 1404 1405 /* Set the datatype */ 1406 #if defined(PETSC_USE_COMPLEX) 1407 lu->datatype = MPI_DOUBLE_COMPLEX; 1408 #else 1409 lu->datatype = MPI_DOUBLE; 1410 #endif 1411 1412 ierr = PLA_Matrix_create(lu->datatype,M,A->cmap->N,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A);CHKERRQ(ierr); 1413 1414 1415 lu->pla_solved = PETSC_FALSE; /* MatSolve_Plapack() is called yet */ 1416 lu->mstruct = DIFFERENT_NONZERO_PATTERN; 1417 1418 if (ftype == MAT_FACTOR_LU) { 1419 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_MPIDense; 1420 } else if (ftype == MAT_FACTOR_CHOLESKY) { 1421 (*F)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MPIDense; 1422 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No incomplete factorizations for dense matrices"); 1423 (*F)->factortype = ftype; 1424 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mpidense_plapack",MatFactorGetSolverPackage_mpidense_plapack);CHKERRQ(ierr); 1425 PetscFunctionReturn(0); 1426 } 1427 EXTERN_C_END 1428 #endif 1429 1430 #undef __FUNCT__ 1431 #define __FUNCT__ "MatGetFactor_mpidense_petsc" 1432 PetscErrorCode MatGetFactor_mpidense_petsc(Mat A,MatFactorType ftype,Mat *F) 1433 { 1434 #if defined(PETSC_HAVE_PLAPACK) 1435 PetscErrorCode ierr; 1436 #endif 1437 1438 PetscFunctionBegin; 1439 #if defined(PETSC_HAVE_PLAPACK) 1440 ierr = MatGetFactor_mpidense_plapack(A,ftype,F);CHKERRQ(ierr); 1441 #else 1442 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix format %s uses PLAPACK direct solver. Install PLAPACK",((PetscObject)A)->type_name); 1443 #endif 1444 PetscFunctionReturn(0); 1445 } 1446 1447 #undef __FUNCT__ 1448 #define __FUNCT__ "MatAXPY_MPIDense" 1449 PetscErrorCode MatAXPY_MPIDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 1450 { 1451 PetscErrorCode ierr; 1452 Mat_MPIDense *A = (Mat_MPIDense*)Y->data, *B = (Mat_MPIDense*)X->data; 1453 1454 PetscFunctionBegin; 1455 ierr = MatAXPY(A->A,alpha,B->A,str);CHKERRQ(ierr); 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #undef __FUNCT__ 1460 #define __FUNCT__ "MatConjugate_MPIDense" 1461 PetscErrorCode MatConjugate_MPIDense(Mat mat) 1462 { 1463 Mat_MPIDense *a = (Mat_MPIDense *)mat->data; 1464 PetscErrorCode ierr; 1465 1466 PetscFunctionBegin; 1467 ierr = MatConjugate(a->A);CHKERRQ(ierr); 1468 PetscFunctionReturn(0); 1469 } 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "MatRealPart_MPIDense" 1473 PetscErrorCode MatRealPart_MPIDense(Mat A) 1474 { 1475 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1476 PetscErrorCode ierr; 1477 1478 PetscFunctionBegin; 1479 ierr = MatRealPart(a->A);CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 #undef __FUNCT__ 1484 #define __FUNCT__ "MatImaginaryPart_MPIDense" 1485 PetscErrorCode MatImaginaryPart_MPIDense(Mat A) 1486 { 1487 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1488 PetscErrorCode ierr; 1489 1490 PetscFunctionBegin; 1491 ierr = MatImaginaryPart(a->A);CHKERRQ(ierr); 1492 PetscFunctionReturn(0); 1493 } 1494 1495 extern PetscErrorCode MatGetColumnNorms_SeqDense(Mat,NormType,PetscReal*); 1496 #undef __FUNCT__ 1497 #define __FUNCT__ "MatGetColumnNorms_MPIDense" 1498 PetscErrorCode MatGetColumnNorms_MPIDense(Mat A,NormType type,PetscReal *norms) 1499 { 1500 PetscErrorCode ierr; 1501 PetscInt i,n; 1502 Mat_MPIDense *a = (Mat_MPIDense*) A->data; 1503 PetscReal *work; 1504 1505 PetscFunctionBegin; 1506 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 1507 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 1508 ierr = MatGetColumnNorms_SeqDense(a->A,type,work);CHKERRQ(ierr); 1509 if (type == NORM_2) { 1510 for (i=0; i<n; i++) work[i] *= work[i]; 1511 } 1512 if (type == NORM_INFINITY) { 1513 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_MAX,A->hdr.comm);CHKERRQ(ierr); 1514 } else { 1515 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_SUM,A->hdr.comm);CHKERRQ(ierr); 1516 } 1517 ierr = PetscFree(work);CHKERRQ(ierr); 1518 if (type == NORM_2) { 1519 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1520 } 1521 PetscFunctionReturn(0); 1522 } 1523 1524 /* -------------------------------------------------------------------*/ 1525 static struct _MatOps MatOps_Values = {MatSetValues_MPIDense, 1526 MatGetRow_MPIDense, 1527 MatRestoreRow_MPIDense, 1528 MatMult_MPIDense, 1529 /* 4*/ MatMultAdd_MPIDense, 1530 MatMultTranspose_MPIDense, 1531 MatMultTransposeAdd_MPIDense, 1532 0, 1533 0, 1534 0, 1535 /*10*/ 0, 1536 0, 1537 0, 1538 0, 1539 MatTranspose_MPIDense, 1540 /*15*/ MatGetInfo_MPIDense, 1541 MatEqual_MPIDense, 1542 MatGetDiagonal_MPIDense, 1543 MatDiagonalScale_MPIDense, 1544 MatNorm_MPIDense, 1545 /*20*/ MatAssemblyBegin_MPIDense, 1546 MatAssemblyEnd_MPIDense, 1547 MatSetOption_MPIDense, 1548 MatZeroEntries_MPIDense, 1549 /*24*/ MatZeroRows_MPIDense, 1550 0, 1551 0, 1552 0, 1553 0, 1554 /*29*/ MatSetUpPreallocation_MPIDense, 1555 0, 1556 0, 1557 MatGetArray_MPIDense, 1558 MatRestoreArray_MPIDense, 1559 /*34*/ MatDuplicate_MPIDense, 1560 0, 1561 0, 1562 0, 1563 0, 1564 /*39*/ MatAXPY_MPIDense, 1565 MatGetSubMatrices_MPIDense, 1566 0, 1567 MatGetValues_MPIDense, 1568 0, 1569 /*44*/ 0, 1570 MatScale_MPIDense, 1571 0, 1572 0, 1573 0, 1574 /*49*/ 0, 1575 0, 1576 0, 1577 0, 1578 0, 1579 /*54*/ 0, 1580 0, 1581 0, 1582 0, 1583 0, 1584 /*59*/ MatGetSubMatrix_MPIDense, 1585 MatDestroy_MPIDense, 1586 MatView_MPIDense, 1587 0, 1588 0, 1589 /*64*/ 0, 1590 0, 1591 0, 1592 0, 1593 0, 1594 /*69*/ 0, 1595 0, 1596 0, 1597 0, 1598 0, 1599 /*74*/ 0, 1600 0, 1601 0, 1602 0, 1603 0, 1604 /*79*/ 0, 1605 0, 1606 0, 1607 0, 1608 /*83*/ MatLoad_MPIDense, 1609 0, 1610 0, 1611 0, 1612 0, 1613 0, 1614 /*89*/ 1615 #if defined(PETSC_HAVE_PLAPACK) 1616 MatMatMult_MPIDense_MPIDense, 1617 MatMatMultSymbolic_MPIDense_MPIDense, 1618 MatMatMultNumeric_MPIDense_MPIDense, 1619 #else 1620 0, 1621 0, 1622 0, 1623 #endif 1624 0, 1625 0, 1626 /*94*/ 0, 1627 0, 1628 0, 1629 0, 1630 0, 1631 /*99*/ 0, 1632 0, 1633 0, 1634 MatConjugate_MPIDense, 1635 0, 1636 /*104*/0, 1637 MatRealPart_MPIDense, 1638 MatImaginaryPart_MPIDense, 1639 0, 1640 0, 1641 /*109*/0, 1642 0, 1643 0, 1644 0, 1645 0, 1646 /*114*/0, 1647 0, 1648 0, 1649 0, 1650 0, 1651 /*119*/0, 1652 0, 1653 0, 1654 0, 1655 0, 1656 /*124*/0, 1657 MatGetColumnNorms_MPIDense 1658 }; 1659 1660 EXTERN_C_BEGIN 1661 #undef __FUNCT__ 1662 #define __FUNCT__ "MatMPIDenseSetPreallocation_MPIDense" 1663 PetscErrorCode MatMPIDenseSetPreallocation_MPIDense(Mat mat,PetscScalar *data) 1664 { 1665 Mat_MPIDense *a; 1666 PetscErrorCode ierr; 1667 1668 PetscFunctionBegin; 1669 mat->preallocated = PETSC_TRUE; 1670 /* Note: For now, when data is specified above, this assumes the user correctly 1671 allocates the local dense storage space. We should add error checking. */ 1672 1673 a = (Mat_MPIDense*)mat->data; 1674 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 1675 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 1676 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 1677 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 1678 a->nvec = mat->cmap->n; 1679 1680 ierr = MatCreate(PETSC_COMM_SELF,&a->A);CHKERRQ(ierr); 1681 ierr = MatSetSizes(a->A,mat->rmap->n,mat->cmap->N,mat->rmap->n,mat->cmap->N);CHKERRQ(ierr); 1682 ierr = MatSetType(a->A,MATSEQDENSE);CHKERRQ(ierr); 1683 ierr = MatSeqDenseSetPreallocation(a->A,data);CHKERRQ(ierr); 1684 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1685 PetscFunctionReturn(0); 1686 } 1687 EXTERN_C_END 1688 1689 /*MC 1690 MATSOLVERPLAPACK = "mpidense" - Parallel LU and Cholesky factorization for MATMPIDENSE matrices 1691 1692 run ./configure with the option --download-plapack 1693 1694 1695 Options Database Keys: 1696 . -mat_plapack_nprows <n> - number of rows in processor partition 1697 . -mat_plapack_npcols <n> - number of columns in processor partition 1698 . -mat_plapack_nb <n> - block size of template vector 1699 . -mat_plapack_nb_alg <n> - algorithmic block size 1700 - -mat_plapack_ckerror <n> - error checking flag 1701 1702 Level: intermediate 1703 1704 .seealso: MatCreateMPIDense(), MATDENSE, MATSEQDENSE, PCFactorSetSolverPackage(), MatSolverPackage 1705 1706 M*/ 1707 1708 EXTERN_C_BEGIN 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "MatCreate_MPIDense" 1711 PetscErrorCode MatCreate_MPIDense(Mat mat) 1712 { 1713 Mat_MPIDense *a; 1714 PetscErrorCode ierr; 1715 1716 PetscFunctionBegin; 1717 ierr = PetscNewLog(mat,Mat_MPIDense,&a);CHKERRQ(ierr); 1718 mat->data = (void*)a; 1719 ierr = PetscMemcpy(mat->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1720 1721 mat->insertmode = NOT_SET_VALUES; 1722 ierr = MPI_Comm_rank(((PetscObject)mat)->comm,&a->rank);CHKERRQ(ierr); 1723 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&a->size);CHKERRQ(ierr); 1724 1725 /* build cache for off array entries formed */ 1726 a->donotstash = PETSC_FALSE; 1727 ierr = MatStashCreate_Private(((PetscObject)mat)->comm,1,&mat->stash);CHKERRQ(ierr); 1728 1729 /* stuff used for matrix vector multiply */ 1730 a->lvec = 0; 1731 a->Mvctx = 0; 1732 a->roworiented = PETSC_TRUE; 1733 1734 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetDiagonalBlock_C", 1735 "MatGetDiagonalBlock_MPIDense", 1736 MatGetDiagonalBlock_MPIDense);CHKERRQ(ierr); 1737 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMPIDenseSetPreallocation_C", 1738 "MatMPIDenseSetPreallocation_MPIDense", 1739 MatMPIDenseSetPreallocation_MPIDense);CHKERRQ(ierr); 1740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_mpiaij_mpidense_C", 1741 "MatMatMult_MPIAIJ_MPIDense", 1742 MatMatMult_MPIAIJ_MPIDense);CHKERRQ(ierr); 1743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_mpiaij_mpidense_C", 1744 "MatMatMultSymbolic_MPIAIJ_MPIDense", 1745 MatMatMultSymbolic_MPIAIJ_MPIDense);CHKERRQ(ierr); 1746 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_mpiaij_mpidense_C", 1747 "MatMatMultNumeric_MPIAIJ_MPIDense", 1748 MatMatMultNumeric_MPIAIJ_MPIDense);CHKERRQ(ierr); 1749 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_petsc_C", 1750 "MatGetFactor_mpidense_petsc", 1751 MatGetFactor_mpidense_petsc);CHKERRQ(ierr); 1752 #if defined(PETSC_HAVE_PLAPACK) 1753 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_plapack_C", 1754 "MatGetFactor_mpidense_plapack", 1755 MatGetFactor_mpidense_plapack);CHKERRQ(ierr); 1756 ierr = PetscPLAPACKInitializePackage(((PetscObject)mat)->comm);CHKERRQ(ierr); 1757 #endif 1758 ierr = PetscObjectChangeTypeName((PetscObject)mat,MATMPIDENSE);CHKERRQ(ierr); 1759 1760 PetscFunctionReturn(0); 1761 } 1762 EXTERN_C_END 1763 1764 /*MC 1765 MATDENSE - MATDENSE = "dense" - A matrix type to be used for dense matrices. 1766 1767 This matrix type is identical to MATSEQDENSE when constructed with a single process communicator, 1768 and MATMPIDENSE otherwise. 1769 1770 Options Database Keys: 1771 . -mat_type dense - sets the matrix type to "dense" during a call to MatSetFromOptions() 1772 1773 Level: beginner 1774 1775 1776 .seealso: MatCreateMPIDense,MATSEQDENSE,MATMPIDENSE 1777 M*/ 1778 1779 #undef __FUNCT__ 1780 #define __FUNCT__ "MatMPIDenseSetPreallocation" 1781 /*@C 1782 MatMPIDenseSetPreallocation - Sets the array used to store the matrix entries 1783 1784 Not collective 1785 1786 Input Parameters: 1787 . A - the matrix 1788 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1789 to control all matrix memory allocation. 1790 1791 Notes: 1792 The dense format is fully compatible with standard Fortran 77 1793 storage by columns. 1794 1795 The data input variable is intended primarily for Fortran programmers 1796 who wish to allocate their own matrix memory space. Most users should 1797 set data=PETSC_NULL. 1798 1799 Level: intermediate 1800 1801 .keywords: matrix,dense, parallel 1802 1803 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1804 @*/ 1805 PetscErrorCode MatMPIDenseSetPreallocation(Mat mat,PetscScalar *data) 1806 { 1807 PetscErrorCode ierr; 1808 1809 PetscFunctionBegin; 1810 ierr = PetscTryMethod(mat,"MatMPIDenseSetPreallocation_C",(Mat,PetscScalar *),(mat,data));CHKERRQ(ierr); 1811 PetscFunctionReturn(0); 1812 } 1813 1814 #undef __FUNCT__ 1815 #define __FUNCT__ "MatCreateMPIDense" 1816 /*@C 1817 MatCreateMPIDense - Creates a sparse parallel matrix in dense format. 1818 1819 Collective on MPI_Comm 1820 1821 Input Parameters: 1822 + comm - MPI communicator 1823 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1824 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 1825 . M - number of global rows (or PETSC_DECIDE to have calculated if m is given) 1826 . N - number of global columns (or PETSC_DECIDE to have calculated if n is given) 1827 - data - optional location of matrix data. Set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users) for PETSc 1828 to control all matrix memory allocation. 1829 1830 Output Parameter: 1831 . A - the matrix 1832 1833 Notes: 1834 The dense format is fully compatible with standard Fortran 77 1835 storage by columns. 1836 1837 The data input variable is intended primarily for Fortran programmers 1838 who wish to allocate their own matrix memory space. Most users should 1839 set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users). 1840 1841 The user MUST specify either the local or global matrix dimensions 1842 (possibly both). 1843 1844 Level: intermediate 1845 1846 .keywords: matrix,dense, parallel 1847 1848 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1849 @*/ 1850 PetscErrorCode MatCreateMPIDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar *data,Mat *A) 1851 { 1852 PetscErrorCode ierr; 1853 PetscMPIInt size; 1854 1855 PetscFunctionBegin; 1856 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1857 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1858 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1859 if (size > 1) { 1860 ierr = MatSetType(*A,MATMPIDENSE);CHKERRQ(ierr); 1861 ierr = MatMPIDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1862 } else { 1863 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1864 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1865 } 1866 PetscFunctionReturn(0); 1867 } 1868 1869 #undef __FUNCT__ 1870 #define __FUNCT__ "MatDuplicate_MPIDense" 1871 static PetscErrorCode MatDuplicate_MPIDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 1872 { 1873 Mat mat; 1874 Mat_MPIDense *a,*oldmat = (Mat_MPIDense*)A->data; 1875 PetscErrorCode ierr; 1876 1877 PetscFunctionBegin; 1878 *newmat = 0; 1879 ierr = MatCreate(((PetscObject)A)->comm,&mat);CHKERRQ(ierr); 1880 ierr = MatSetSizes(mat,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1881 ierr = MatSetType(mat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1882 a = (Mat_MPIDense*)mat->data; 1883 ierr = PetscMemcpy(mat->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1884 1885 mat->factortype = A->factortype; 1886 mat->assembled = PETSC_TRUE; 1887 mat->preallocated = PETSC_TRUE; 1888 1889 a->size = oldmat->size; 1890 a->rank = oldmat->rank; 1891 mat->insertmode = NOT_SET_VALUES; 1892 a->nvec = oldmat->nvec; 1893 a->donotstash = oldmat->donotstash; 1894 1895 ierr = PetscLayoutReference(A->rmap,&mat->rmap);CHKERRQ(ierr); 1896 ierr = PetscLayoutReference(A->cmap,&mat->cmap);CHKERRQ(ierr); 1897 1898 ierr = MatSetUpMultiply_MPIDense(mat);CHKERRQ(ierr); 1899 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 1900 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1901 1902 *newmat = mat; 1903 PetscFunctionReturn(0); 1904 } 1905 1906 #undef __FUNCT__ 1907 #define __FUNCT__ "MatLoad_MPIDense_DenseInFile" 1908 PetscErrorCode MatLoad_MPIDense_DenseInFile(MPI_Comm comm,PetscInt fd,PetscInt M,PetscInt N,Mat newmat,PetscInt sizesset) 1909 { 1910 PetscErrorCode ierr; 1911 PetscMPIInt rank,size; 1912 PetscInt *rowners,i,m,nz,j; 1913 PetscScalar *array,*vals,*vals_ptr; 1914 1915 PetscFunctionBegin; 1916 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1917 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1918 1919 /* determine ownership of all rows */ 1920 if (newmat->rmap->n < 0) m = M/size + ((M % size) > rank); 1921 else m = newmat->rmap->n; 1922 ierr = PetscMalloc((size+2)*sizeof(PetscInt),&rowners);CHKERRQ(ierr); 1923 ierr = MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 1924 rowners[0] = 0; 1925 for (i=2; i<=size; i++) { 1926 rowners[i] += rowners[i-1]; 1927 } 1928 1929 if (!sizesset) { 1930 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 1931 } 1932 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1933 ierr = MatGetArray(newmat,&array);CHKERRQ(ierr); 1934 1935 if (!rank) { 1936 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1937 1938 /* read in my part of the matrix numerical values */ 1939 ierr = PetscBinaryRead(fd,vals,m*N,PETSC_SCALAR);CHKERRQ(ierr); 1940 1941 /* insert into matrix-by row (this is why cannot directly read into array */ 1942 vals_ptr = vals; 1943 for (i=0; i<m; i++) { 1944 for (j=0; j<N; j++) { 1945 array[i + j*m] = *vals_ptr++; 1946 } 1947 } 1948 1949 /* read in other processors and ship out */ 1950 for (i=1; i<size; i++) { 1951 nz = (rowners[i+1] - rowners[i])*N; 1952 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 1953 ierr = MPILong_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1954 } 1955 } else { 1956 /* receive numeric values */ 1957 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1958 1959 /* receive message of values*/ 1960 ierr = MPILong_Recv(vals,m*N,MPIU_SCALAR,0,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1961 1962 /* insert into matrix-by row (this is why cannot directly read into array */ 1963 vals_ptr = vals; 1964 for (i=0; i<m; i++) { 1965 for (j=0; j<N; j++) { 1966 array[i + j*m] = *vals_ptr++; 1967 } 1968 } 1969 } 1970 ierr = PetscFree(rowners);CHKERRQ(ierr); 1971 ierr = PetscFree(vals);CHKERRQ(ierr); 1972 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1973 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1974 PetscFunctionReturn(0); 1975 } 1976 1977 #undef __FUNCT__ 1978 #define __FUNCT__ "MatLoad_MPIDense" 1979 PetscErrorCode MatLoad_MPIDense(Mat newmat,PetscViewer viewer) 1980 { 1981 PetscScalar *vals,*svals; 1982 MPI_Comm comm = ((PetscObject)viewer)->comm; 1983 MPI_Status status; 1984 PetscMPIInt rank,size,tag = ((PetscObject)viewer)->tag,*rowners,*sndcounts,m,maxnz; 1985 PetscInt header[4],*rowlengths = 0,M,N,*cols; 1986 PetscInt *ourlens,*procsnz = 0,*offlens,jj,*mycols,*smycols; 1987 PetscInt i,nz,j,rstart,rend,sizesset=1,grows,gcols; 1988 int fd; 1989 PetscErrorCode ierr; 1990 1991 PetscFunctionBegin; 1992 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1993 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1994 if (!rank) { 1995 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1996 ierr = PetscBinaryRead(fd,(char *)header,4,PETSC_INT);CHKERRQ(ierr); 1997 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object"); 1998 } 1999 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) sizesset = 0; 2000 2001 ierr = MPI_Bcast(header+1,3,MPIU_INT,0,comm);CHKERRQ(ierr); 2002 M = header[1]; N = header[2]; nz = header[3]; 2003 2004 /* If global rows/cols are set to PETSC_DECIDE, set it to the sizes given in the file */ 2005 if (sizesset && newmat->rmap->N < 0) newmat->rmap->N = M; 2006 if (sizesset && newmat->cmap->N < 0) newmat->cmap->N = N; 2007 2008 /* If global sizes are set, check if they are consistent with that given in the file */ 2009 if (sizesset) { 2010 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 2011 } 2012 if (sizesset && newmat->rmap->N != grows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Inconsistent # of rows:Matrix in file has (%d) and input matrix has (%d)",M,grows); 2013 if (sizesset && newmat->cmap->N != gcols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Inconsistent # of cols:Matrix in file has (%d) and input matrix has (%d)",N,gcols); 2014 2015 /* 2016 Handle case where matrix is stored on disk as a dense matrix 2017 */ 2018 if (nz == MATRIX_BINARY_FORMAT_DENSE) { 2019 ierr = MatLoad_MPIDense_DenseInFile(comm,fd,M,N,newmat,sizesset);CHKERRQ(ierr); 2020 PetscFunctionReturn(0); 2021 } 2022 2023 /* determine ownership of all rows */ 2024 if (newmat->rmap->n < 0) m = PetscMPIIntCast(M/size + ((M % size) > rank)); 2025 else m = PetscMPIIntCast(newmat->rmap->n); 2026 ierr = PetscMalloc((size+2)*sizeof(PetscMPIInt),&rowners);CHKERRQ(ierr); 2027 ierr = MPI_Allgather(&m,1,MPI_INT,rowners+1,1,MPI_INT,comm);CHKERRQ(ierr); 2028 rowners[0] = 0; 2029 for (i=2; i<=size; i++) { 2030 rowners[i] += rowners[i-1]; 2031 } 2032 rstart = rowners[rank]; 2033 rend = rowners[rank+1]; 2034 2035 /* distribute row lengths to all processors */ 2036 ierr = PetscMalloc2(rend-rstart,PetscInt,&ourlens,rend-rstart,PetscInt,&offlens);CHKERRQ(ierr); 2037 if (!rank) { 2038 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2039 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2040 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&sndcounts);CHKERRQ(ierr); 2041 for (i=0; i<size; i++) sndcounts[i] = rowners[i+1] - rowners[i]; 2042 ierr = MPI_Scatterv(rowlengths,sndcounts,rowners,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2043 ierr = PetscFree(sndcounts);CHKERRQ(ierr); 2044 } else { 2045 ierr = MPI_Scatterv(0,0,0,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2046 } 2047 2048 if (!rank) { 2049 /* calculate the number of nonzeros on each processor */ 2050 ierr = PetscMalloc(size*sizeof(PetscInt),&procsnz);CHKERRQ(ierr); 2051 ierr = PetscMemzero(procsnz,size*sizeof(PetscInt));CHKERRQ(ierr); 2052 for (i=0; i<size; i++) { 2053 for (j=rowners[i]; j< rowners[i+1]; j++) { 2054 procsnz[i] += rowlengths[j]; 2055 } 2056 } 2057 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2058 2059 /* determine max buffer needed and allocate it */ 2060 maxnz = 0; 2061 for (i=0; i<size; i++) { 2062 maxnz = PetscMax(maxnz,procsnz[i]); 2063 } 2064 ierr = PetscMalloc(maxnz*sizeof(PetscInt),&cols);CHKERRQ(ierr); 2065 2066 /* read in my part of the matrix column indices */ 2067 nz = procsnz[0]; 2068 ierr = PetscMalloc(nz*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2069 ierr = PetscBinaryRead(fd,mycols,nz,PETSC_INT);CHKERRQ(ierr); 2070 2071 /* read in every one elses and ship off */ 2072 for (i=1; i<size; i++) { 2073 nz = procsnz[i]; 2074 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 2075 ierr = MPI_Send(cols,nz,MPIU_INT,i,tag,comm);CHKERRQ(ierr); 2076 } 2077 ierr = PetscFree(cols);CHKERRQ(ierr); 2078 } else { 2079 /* determine buffer space needed for message */ 2080 nz = 0; 2081 for (i=0; i<m; i++) { 2082 nz += ourlens[i]; 2083 } 2084 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2085 2086 /* receive message of column indices*/ 2087 ierr = MPI_Recv(mycols,nz,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 2088 ierr = MPI_Get_count(&status,MPIU_INT,&maxnz);CHKERRQ(ierr); 2089 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2090 } 2091 2092 /* loop over local rows, determining number of off diagonal entries */ 2093 ierr = PetscMemzero(offlens,m*sizeof(PetscInt));CHKERRQ(ierr); 2094 jj = 0; 2095 for (i=0; i<m; i++) { 2096 for (j=0; j<ourlens[i]; j++) { 2097 if (mycols[jj] < rstart || mycols[jj] >= rend) offlens[i]++; 2098 jj++; 2099 } 2100 } 2101 2102 /* create our matrix */ 2103 for (i=0; i<m; i++) { 2104 ourlens[i] -= offlens[i]; 2105 } 2106 2107 if (!sizesset) { 2108 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2109 } 2110 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 2111 for (i=0; i<m; i++) { 2112 ourlens[i] += offlens[i]; 2113 } 2114 2115 if (!rank) { 2116 ierr = PetscMalloc(maxnz*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2117 2118 /* read in my part of the matrix numerical values */ 2119 nz = procsnz[0]; 2120 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2121 2122 /* insert into matrix */ 2123 jj = rstart; 2124 smycols = mycols; 2125 svals = vals; 2126 for (i=0; i<m; i++) { 2127 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2128 smycols += ourlens[i]; 2129 svals += ourlens[i]; 2130 jj++; 2131 } 2132 2133 /* read in other processors and ship out */ 2134 for (i=1; i<size; i++) { 2135 nz = procsnz[i]; 2136 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2137 ierr = MPI_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 2138 } 2139 ierr = PetscFree(procsnz);CHKERRQ(ierr); 2140 } else { 2141 /* receive numeric values */ 2142 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2143 2144 /* receive message of values*/ 2145 ierr = MPI_Recv(vals,nz,MPIU_SCALAR,0,((PetscObject)newmat)->tag,comm,&status);CHKERRQ(ierr); 2146 ierr = MPI_Get_count(&status,MPIU_SCALAR,&maxnz);CHKERRQ(ierr); 2147 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2148 2149 /* insert into matrix */ 2150 jj = rstart; 2151 smycols = mycols; 2152 svals = vals; 2153 for (i=0; i<m; i++) { 2154 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2155 smycols += ourlens[i]; 2156 svals += ourlens[i]; 2157 jj++; 2158 } 2159 } 2160 ierr = PetscFree2(ourlens,offlens);CHKERRQ(ierr); 2161 ierr = PetscFree(vals);CHKERRQ(ierr); 2162 ierr = PetscFree(mycols);CHKERRQ(ierr); 2163 ierr = PetscFree(rowners);CHKERRQ(ierr); 2164 2165 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2166 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2167 PetscFunctionReturn(0); 2168 } 2169 2170 #undef __FUNCT__ 2171 #define __FUNCT__ "MatEqual_MPIDense" 2172 PetscErrorCode MatEqual_MPIDense(Mat A,Mat B,PetscBool *flag) 2173 { 2174 Mat_MPIDense *matB = (Mat_MPIDense*)B->data,*matA = (Mat_MPIDense*)A->data; 2175 Mat a,b; 2176 PetscBool flg; 2177 PetscErrorCode ierr; 2178 2179 PetscFunctionBegin; 2180 a = matA->A; 2181 b = matB->A; 2182 ierr = MatEqual(a,b,&flg);CHKERRQ(ierr); 2183 ierr = MPI_Allreduce(&flg,flag,1,MPI_INT,MPI_LAND,((PetscObject)A)->comm);CHKERRQ(ierr); 2184 PetscFunctionReturn(0); 2185 } 2186 2187 #if defined(PETSC_HAVE_PLAPACK) 2188 2189 #undef __FUNCT__ 2190 #define __FUNCT__ "PetscPLAPACKFinalizePackage" 2191 /*@C 2192 PetscPLAPACKFinalizePackage - This function destroys everything in the Petsc interface to PLAPACK. 2193 Level: developer 2194 2195 .keywords: Petsc, destroy, package, PLAPACK 2196 .seealso: PetscFinalize() 2197 @*/ 2198 PetscErrorCode PetscPLAPACKFinalizePackage(void) 2199 { 2200 PetscErrorCode ierr; 2201 2202 PetscFunctionBegin; 2203 ierr = PLA_Finalize();CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 #undef __FUNCT__ 2208 #define __FUNCT__ "PetscPLAPACKInitializePackage" 2209 /*@C 2210 PetscPLAPACKInitializePackage - This function initializes everything in the Petsc interface to PLAPACK. It is 2211 called from MatCreate_MPIDense() the first time an MPI dense matrix is called. 2212 2213 Input Parameter: 2214 . comm - the communicator the matrix lives on 2215 2216 Level: developer 2217 2218 Notes: PLAPACK does not have a good fit with MPI communicators; all (parallel) PLAPACK objects have to live in the 2219 same communicator (because there is some global state that is initialized and used for all matrices). In addition if 2220 PLAPACK is initialized (that is the initial matrices created) are on subcommunicators of MPI_COMM_WORLD, these subcommunicators 2221 cannot overlap. 2222 2223 .keywords: Petsc, initialize, package, PLAPACK 2224 .seealso: PetscSysInitializePackage(), PetscInitialize() 2225 @*/ 2226 PetscErrorCode PetscPLAPACKInitializePackage(MPI_Comm comm) 2227 { 2228 PetscMPIInt size; 2229 PetscErrorCode ierr; 2230 2231 PetscFunctionBegin; 2232 if (!PLA_Initialized(PETSC_NULL)) { 2233 2234 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2235 Plapack_nprows = 1; 2236 Plapack_npcols = size; 2237 2238 ierr = PetscOptionsBegin(comm,PETSC_NULL,"PLAPACK Options","Mat");CHKERRQ(ierr); 2239 ierr = PetscOptionsInt("-mat_plapack_nprows","row dimension of 2D processor mesh","None",Plapack_nprows,&Plapack_nprows,PETSC_NULL);CHKERRQ(ierr); 2240 ierr = PetscOptionsInt("-mat_plapack_npcols","column dimension of 2D processor mesh","None",Plapack_npcols,&Plapack_npcols,PETSC_NULL);CHKERRQ(ierr); 2241 #if defined(PETSC_USE_DEBUG) 2242 Plapack_ierror = 3; 2243 #else 2244 Plapack_ierror = 0; 2245 #endif 2246 ierr = PetscOptionsInt("-mat_plapack_ckerror","error checking flag","None",Plapack_ierror,&Plapack_ierror,PETSC_NULL);CHKERRQ(ierr); 2247 if (Plapack_ierror){ 2248 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );CHKERRQ(ierr); 2249 } else { 2250 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );CHKERRQ(ierr); 2251 } 2252 2253 Plapack_nb_alg = 0; 2254 ierr = PetscOptionsInt("-mat_plapack_nb_alg","algorithmic block size","None",Plapack_nb_alg,&Plapack_nb_alg,PETSC_NULL);CHKERRQ(ierr); 2255 if (Plapack_nb_alg) { 2256 ierr = pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,Plapack_nb_alg);CHKERRQ(ierr); 2257 } 2258 PetscOptionsEnd(); 2259 2260 ierr = PLA_Comm_1D_to_2D(comm,Plapack_nprows,Plapack_npcols,&Plapack_comm_2d);CHKERRQ(ierr); 2261 ierr = PLA_Init(Plapack_comm_2d);CHKERRQ(ierr); 2262 ierr = PetscRegisterFinalize(PetscPLAPACKFinalizePackage);CHKERRQ(ierr); 2263 } 2264 PetscFunctionReturn(0); 2265 } 2266 2267 #endif 2268