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