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