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 = PetscObjectTypeCompare((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 = MPIULong_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 = MPIULong_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 = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 701 ierr = PetscObjectTypeCompare((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 = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 794 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 795 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSOCKET,&issocket);CHKERRQ(ierr); 796 ierr = PetscObjectTypeCompare((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,1,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 F->preallocated = PETSC_TRUE; 1349 PetscFunctionReturn(0); 1350 } 1351 1352 /* Note the Petsc r and c permutations are ignored */ 1353 #undef __FUNCT__ 1354 #define __FUNCT__ "MatLUFactorSymbolic_MPIDense" 1355 PetscErrorCode MatLUFactorSymbolic_MPIDense(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1356 { 1357 PetscErrorCode ierr; 1358 PetscInt M = A->rmap->N; 1359 Mat_Plapack *lu; 1360 1361 PetscFunctionBegin; 1362 lu = (Mat_Plapack*)F->spptr; 1363 ierr = PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots);CHKERRQ(ierr); 1364 F->ops->lufactornumeric = MatLUFactorNumeric_MPIDense; 1365 F->preallocated = PETSC_TRUE; 1366 PetscFunctionReturn(0); 1367 } 1368 1369 EXTERN_C_BEGIN 1370 #undef __FUNCT__ 1371 #define __FUNCT__ "MatFactorGetSolverPackage_mpidense_plapack" 1372 PetscErrorCode MatFactorGetSolverPackage_mpidense_plapack(Mat A,const MatSolverPackage *type) 1373 { 1374 PetscFunctionBegin; 1375 *type = MATSOLVERPLAPACK; 1376 PetscFunctionReturn(0); 1377 } 1378 EXTERN_C_END 1379 1380 EXTERN_C_BEGIN 1381 #undef __FUNCT__ 1382 #define __FUNCT__ "MatGetFactor_mpidense_plapack" 1383 PetscErrorCode MatGetFactor_mpidense_plapack(Mat A,MatFactorType ftype,Mat *F) 1384 { 1385 PetscErrorCode ierr; 1386 Mat_Plapack *lu; 1387 PetscMPIInt size; 1388 PetscInt M=A->rmap->N; 1389 1390 PetscFunctionBegin; 1391 /* Create the factorization matrix */ 1392 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 1393 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1394 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 1395 ierr = MatSetUp(*F);CHKERRQ(ierr); 1396 ierr = PetscNewLog(*F,Mat_Plapack,&lu);CHKERRQ(ierr); 1397 (*F)->spptr = (void*)lu; 1398 1399 /* Set default Plapack parameters */ 1400 ierr = MPI_Comm_size(((PetscObject)A)->comm,&size);CHKERRQ(ierr); 1401 lu->nb = M/size; 1402 if (M - lu->nb*size) lu->nb++; /* without cyclic distribution */ 1403 1404 /* Set runtime options */ 1405 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"PLAPACK Options","Mat");CHKERRQ(ierr); 1406 ierr = PetscOptionsInt("-mat_plapack_nb","block size of template vector","None",lu->nb,&lu->nb,PETSC_NULL);CHKERRQ(ierr); 1407 PetscOptionsEnd(); 1408 1409 /* Create object distribution template */ 1410 lu->templ = NULL; 1411 ierr = PLA_Temp_create(lu->nb, 0, &lu->templ);CHKERRQ(ierr); 1412 1413 /* Set the datatype */ 1414 #if defined(PETSC_USE_COMPLEX) 1415 lu->datatype = MPI_DOUBLE_COMPLEX; 1416 #else 1417 lu->datatype = MPI_DOUBLE; 1418 #endif 1419 1420 ierr = PLA_Matrix_create(lu->datatype,M,A->cmap->N,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A);CHKERRQ(ierr); 1421 1422 1423 lu->pla_solved = PETSC_FALSE; /* MatSolve_Plapack() is called yet */ 1424 lu->mstruct = DIFFERENT_NONZERO_PATTERN; 1425 1426 if (ftype == MAT_FACTOR_LU) { 1427 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_MPIDense; 1428 } else if (ftype == MAT_FACTOR_CHOLESKY) { 1429 (*F)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MPIDense; 1430 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No incomplete factorizations for dense matrices"); 1431 (*F)->factortype = ftype; 1432 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mpidense_plapack",MatFactorGetSolverPackage_mpidense_plapack);CHKERRQ(ierr); 1433 PetscFunctionReturn(0); 1434 } 1435 EXTERN_C_END 1436 #endif 1437 1438 #undef __FUNCT__ 1439 #define __FUNCT__ "MatGetFactor_mpidense_petsc" 1440 PetscErrorCode MatGetFactor_mpidense_petsc(Mat A,MatFactorType ftype,Mat *F) 1441 { 1442 #if defined(PETSC_HAVE_PLAPACK) 1443 PetscErrorCode ierr; 1444 #endif 1445 1446 PetscFunctionBegin; 1447 #if defined(PETSC_HAVE_PLAPACK) 1448 ierr = MatGetFactor_mpidense_plapack(A,ftype,F);CHKERRQ(ierr); 1449 #else 1450 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix format %s uses PLAPACK direct solver. Install PLAPACK",((PetscObject)A)->type_name); 1451 #endif 1452 PetscFunctionReturn(0); 1453 } 1454 1455 #undef __FUNCT__ 1456 #define __FUNCT__ "MatAXPY_MPIDense" 1457 PetscErrorCode MatAXPY_MPIDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 1458 { 1459 PetscErrorCode ierr; 1460 Mat_MPIDense *A = (Mat_MPIDense*)Y->data, *B = (Mat_MPIDense*)X->data; 1461 1462 PetscFunctionBegin; 1463 ierr = MatAXPY(A->A,alpha,B->A,str);CHKERRQ(ierr); 1464 PetscFunctionReturn(0); 1465 } 1466 1467 #undef __FUNCT__ 1468 #define __FUNCT__ "MatConjugate_MPIDense" 1469 PetscErrorCode MatConjugate_MPIDense(Mat mat) 1470 { 1471 Mat_MPIDense *a = (Mat_MPIDense *)mat->data; 1472 PetscErrorCode ierr; 1473 1474 PetscFunctionBegin; 1475 ierr = MatConjugate(a->A);CHKERRQ(ierr); 1476 PetscFunctionReturn(0); 1477 } 1478 1479 #undef __FUNCT__ 1480 #define __FUNCT__ "MatRealPart_MPIDense" 1481 PetscErrorCode MatRealPart_MPIDense(Mat A) 1482 { 1483 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1484 PetscErrorCode ierr; 1485 1486 PetscFunctionBegin; 1487 ierr = MatRealPart(a->A);CHKERRQ(ierr); 1488 PetscFunctionReturn(0); 1489 } 1490 1491 #undef __FUNCT__ 1492 #define __FUNCT__ "MatImaginaryPart_MPIDense" 1493 PetscErrorCode MatImaginaryPart_MPIDense(Mat A) 1494 { 1495 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1496 PetscErrorCode ierr; 1497 1498 PetscFunctionBegin; 1499 ierr = MatImaginaryPart(a->A);CHKERRQ(ierr); 1500 PetscFunctionReturn(0); 1501 } 1502 1503 extern PetscErrorCode MatGetColumnNorms_SeqDense(Mat,NormType,PetscReal*); 1504 #undef __FUNCT__ 1505 #define __FUNCT__ "MatGetColumnNorms_MPIDense" 1506 PetscErrorCode MatGetColumnNorms_MPIDense(Mat A,NormType type,PetscReal *norms) 1507 { 1508 PetscErrorCode ierr; 1509 PetscInt i,n; 1510 Mat_MPIDense *a = (Mat_MPIDense*) A->data; 1511 PetscReal *work; 1512 1513 PetscFunctionBegin; 1514 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 1515 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 1516 ierr = MatGetColumnNorms_SeqDense(a->A,type,work);CHKERRQ(ierr); 1517 if (type == NORM_2) { 1518 for (i=0; i<n; i++) work[i] *= work[i]; 1519 } 1520 if (type == NORM_INFINITY) { 1521 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_MAX,A->hdr.comm);CHKERRQ(ierr); 1522 } else { 1523 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_SUM,A->hdr.comm);CHKERRQ(ierr); 1524 } 1525 ierr = PetscFree(work);CHKERRQ(ierr); 1526 if (type == NORM_2) { 1527 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1528 } 1529 PetscFunctionReturn(0); 1530 } 1531 1532 /* -------------------------------------------------------------------*/ 1533 static struct _MatOps MatOps_Values = {MatSetValues_MPIDense, 1534 MatGetRow_MPIDense, 1535 MatRestoreRow_MPIDense, 1536 MatMult_MPIDense, 1537 /* 4*/ MatMultAdd_MPIDense, 1538 MatMultTranspose_MPIDense, 1539 MatMultTransposeAdd_MPIDense, 1540 0, 1541 0, 1542 0, 1543 /*10*/ 0, 1544 0, 1545 0, 1546 0, 1547 MatTranspose_MPIDense, 1548 /*15*/ MatGetInfo_MPIDense, 1549 MatEqual_MPIDense, 1550 MatGetDiagonal_MPIDense, 1551 MatDiagonalScale_MPIDense, 1552 MatNorm_MPIDense, 1553 /*20*/ MatAssemblyBegin_MPIDense, 1554 MatAssemblyEnd_MPIDense, 1555 MatSetOption_MPIDense, 1556 MatZeroEntries_MPIDense, 1557 /*24*/ MatZeroRows_MPIDense, 1558 0, 1559 0, 1560 0, 1561 0, 1562 /*29*/ MatSetUp_MPIDense, 1563 0, 1564 0, 1565 MatGetArray_MPIDense, 1566 MatRestoreArray_MPIDense, 1567 /*34*/ MatDuplicate_MPIDense, 1568 0, 1569 0, 1570 0, 1571 0, 1572 /*39*/ MatAXPY_MPIDense, 1573 MatGetSubMatrices_MPIDense, 1574 0, 1575 MatGetValues_MPIDense, 1576 0, 1577 /*44*/ 0, 1578 MatScale_MPIDense, 1579 0, 1580 0, 1581 0, 1582 /*49*/ 0, 1583 0, 1584 0, 1585 0, 1586 0, 1587 /*54*/ 0, 1588 0, 1589 0, 1590 0, 1591 0, 1592 /*59*/ MatGetSubMatrix_MPIDense, 1593 MatDestroy_MPIDense, 1594 MatView_MPIDense, 1595 0, 1596 0, 1597 /*64*/ 0, 1598 0, 1599 0, 1600 0, 1601 0, 1602 /*69*/ 0, 1603 0, 1604 0, 1605 0, 1606 0, 1607 /*74*/ 0, 1608 0, 1609 0, 1610 0, 1611 0, 1612 /*79*/ 0, 1613 0, 1614 0, 1615 0, 1616 /*83*/ MatLoad_MPIDense, 1617 0, 1618 0, 1619 0, 1620 0, 1621 0, 1622 /*89*/ 1623 #if defined(PETSC_HAVE_PLAPACK) 1624 MatMatMult_MPIDense_MPIDense, 1625 MatMatMultSymbolic_MPIDense_MPIDense, 1626 MatMatMultNumeric_MPIDense_MPIDense, 1627 #else 1628 0, 1629 0, 1630 0, 1631 #endif 1632 0, 1633 0, 1634 /*94*/ 0, 1635 0, 1636 0, 1637 0, 1638 0, 1639 /*99*/ 0, 1640 0, 1641 0, 1642 MatConjugate_MPIDense, 1643 0, 1644 /*104*/0, 1645 MatRealPart_MPIDense, 1646 MatImaginaryPart_MPIDense, 1647 0, 1648 0, 1649 /*109*/0, 1650 0, 1651 0, 1652 0, 1653 0, 1654 /*114*/0, 1655 0, 1656 0, 1657 0, 1658 0, 1659 /*119*/0, 1660 0, 1661 0, 1662 0, 1663 0, 1664 /*124*/0, 1665 MatGetColumnNorms_MPIDense 1666 }; 1667 1668 EXTERN_C_BEGIN 1669 #undef __FUNCT__ 1670 #define __FUNCT__ "MatMPIDenseSetPreallocation_MPIDense" 1671 PetscErrorCode MatMPIDenseSetPreallocation_MPIDense(Mat mat,PetscScalar *data) 1672 { 1673 Mat_MPIDense *a; 1674 PetscErrorCode ierr; 1675 1676 PetscFunctionBegin; 1677 mat->preallocated = PETSC_TRUE; 1678 /* Note: For now, when data is specified above, this assumes the user correctly 1679 allocates the local dense storage space. We should add error checking. */ 1680 1681 a = (Mat_MPIDense*)mat->data; 1682 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 1683 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 1684 a->nvec = mat->cmap->n; 1685 1686 ierr = MatCreate(PETSC_COMM_SELF,&a->A);CHKERRQ(ierr); 1687 ierr = MatSetSizes(a->A,mat->rmap->n,mat->cmap->N,mat->rmap->n,mat->cmap->N);CHKERRQ(ierr); 1688 ierr = MatSetType(a->A,MATSEQDENSE);CHKERRQ(ierr); 1689 ierr = MatSeqDenseSetPreallocation(a->A,data);CHKERRQ(ierr); 1690 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1691 PetscFunctionReturn(0); 1692 } 1693 EXTERN_C_END 1694 1695 /*MC 1696 MATSOLVERPLAPACK = "mpidense" - Parallel LU and Cholesky factorization for MATMPIDENSE matrices 1697 1698 run ./configure with the option --download-plapack 1699 1700 1701 Options Database Keys: 1702 . -mat_plapack_nprows <n> - number of rows in processor partition 1703 . -mat_plapack_npcols <n> - number of columns in processor partition 1704 . -mat_plapack_nb <n> - block size of template vector 1705 . -mat_plapack_nb_alg <n> - algorithmic block size 1706 - -mat_plapack_ckerror <n> - error checking flag 1707 1708 Level: intermediate 1709 1710 .seealso: MatCreateDense(), MATDENSE, MATSEQDENSE, PCFactorSetSolverPackage(), MatSolverPackage 1711 1712 M*/ 1713 1714 EXTERN_C_BEGIN 1715 #undef __FUNCT__ 1716 #define __FUNCT__ "MatCreate_MPIDense" 1717 PetscErrorCode MatCreate_MPIDense(Mat mat) 1718 { 1719 Mat_MPIDense *a; 1720 PetscErrorCode ierr; 1721 1722 PetscFunctionBegin; 1723 ierr = PetscNewLog(mat,Mat_MPIDense,&a);CHKERRQ(ierr); 1724 mat->data = (void*)a; 1725 ierr = PetscMemcpy(mat->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1726 1727 mat->insertmode = NOT_SET_VALUES; 1728 ierr = MPI_Comm_rank(((PetscObject)mat)->comm,&a->rank);CHKERRQ(ierr); 1729 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&a->size);CHKERRQ(ierr); 1730 1731 /* build cache for off array entries formed */ 1732 a->donotstash = PETSC_FALSE; 1733 ierr = MatStashCreate_Private(((PetscObject)mat)->comm,1,&mat->stash);CHKERRQ(ierr); 1734 1735 /* stuff used for matrix vector multiply */ 1736 a->lvec = 0; 1737 a->Mvctx = 0; 1738 a->roworiented = PETSC_TRUE; 1739 1740 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetDiagonalBlock_C", 1741 "MatGetDiagonalBlock_MPIDense", 1742 MatGetDiagonalBlock_MPIDense);CHKERRQ(ierr); 1743 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMPIDenseSetPreallocation_C", 1744 "MatMPIDenseSetPreallocation_MPIDense", 1745 MatMPIDenseSetPreallocation_MPIDense);CHKERRQ(ierr); 1746 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_mpiaij_mpidense_C", 1747 "MatMatMult_MPIAIJ_MPIDense", 1748 MatMatMult_MPIAIJ_MPIDense);CHKERRQ(ierr); 1749 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_mpiaij_mpidense_C", 1750 "MatMatMultSymbolic_MPIAIJ_MPIDense", 1751 MatMatMultSymbolic_MPIAIJ_MPIDense);CHKERRQ(ierr); 1752 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_mpiaij_mpidense_C", 1753 "MatMatMultNumeric_MPIAIJ_MPIDense", 1754 MatMatMultNumeric_MPIAIJ_MPIDense);CHKERRQ(ierr); 1755 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_petsc_C", 1756 "MatGetFactor_mpidense_petsc", 1757 MatGetFactor_mpidense_petsc);CHKERRQ(ierr); 1758 #if defined(PETSC_HAVE_PLAPACK) 1759 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_plapack_C", 1760 "MatGetFactor_mpidense_plapack", 1761 MatGetFactor_mpidense_plapack);CHKERRQ(ierr); 1762 ierr = PetscPLAPACKInitializePackage(((PetscObject)mat)->comm);CHKERRQ(ierr); 1763 #endif 1764 ierr = PetscObjectChangeTypeName((PetscObject)mat,MATMPIDENSE);CHKERRQ(ierr); 1765 1766 PetscFunctionReturn(0); 1767 } 1768 EXTERN_C_END 1769 1770 /*MC 1771 MATDENSE - MATDENSE = "dense" - A matrix type to be used for dense matrices. 1772 1773 This matrix type is identical to MATSEQDENSE when constructed with a single process communicator, 1774 and MATMPIDENSE otherwise. 1775 1776 Options Database Keys: 1777 . -mat_type dense - sets the matrix type to "dense" during a call to MatSetFromOptions() 1778 1779 Level: beginner 1780 1781 1782 .seealso: MatCreateMPIDense,MATSEQDENSE,MATMPIDENSE 1783 M*/ 1784 1785 #undef __FUNCT__ 1786 #define __FUNCT__ "MatMPIDenseSetPreallocation" 1787 /*@C 1788 MatMPIDenseSetPreallocation - Sets the array used to store the matrix entries 1789 1790 Not collective 1791 1792 Input Parameters: 1793 . A - the matrix 1794 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1795 to control all matrix memory allocation. 1796 1797 Notes: 1798 The dense format is fully compatible with standard Fortran 77 1799 storage by columns. 1800 1801 The data input variable is intended primarily for Fortran programmers 1802 who wish to allocate their own matrix memory space. Most users should 1803 set data=PETSC_NULL. 1804 1805 Level: intermediate 1806 1807 .keywords: matrix,dense, parallel 1808 1809 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1810 @*/ 1811 PetscErrorCode MatMPIDenseSetPreallocation(Mat mat,PetscScalar *data) 1812 { 1813 PetscErrorCode ierr; 1814 1815 PetscFunctionBegin; 1816 ierr = PetscTryMethod(mat,"MatMPIDenseSetPreallocation_C",(Mat,PetscScalar *),(mat,data));CHKERRQ(ierr); 1817 PetscFunctionReturn(0); 1818 } 1819 1820 #undef __FUNCT__ 1821 #define __FUNCT__ "MatCreateDense" 1822 /*@C 1823 MatCreateDense - Creates a parallel matrix in dense format. 1824 1825 Collective on MPI_Comm 1826 1827 Input Parameters: 1828 + comm - MPI communicator 1829 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1830 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 1831 . M - number of global rows (or PETSC_DECIDE to have calculated if m is given) 1832 . N - number of global columns (or PETSC_DECIDE to have calculated if n is given) 1833 - data - optional location of matrix data. Set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users) for PETSc 1834 to control all matrix memory allocation. 1835 1836 Output Parameter: 1837 . A - the matrix 1838 1839 Notes: 1840 The dense format is fully compatible with standard Fortran 77 1841 storage by columns. 1842 1843 The data input variable is intended primarily for Fortran programmers 1844 who wish to allocate their own matrix memory space. Most users should 1845 set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users). 1846 1847 The user MUST specify either the local or global matrix dimensions 1848 (possibly both). 1849 1850 Level: intermediate 1851 1852 .keywords: matrix,dense, parallel 1853 1854 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1855 @*/ 1856 PetscErrorCode MatCreateDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar *data,Mat *A) 1857 { 1858 PetscErrorCode ierr; 1859 PetscMPIInt size; 1860 1861 PetscFunctionBegin; 1862 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1863 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1864 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1865 if (size > 1) { 1866 ierr = MatSetType(*A,MATMPIDENSE);CHKERRQ(ierr); 1867 ierr = MatMPIDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1868 } else { 1869 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1870 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1871 } 1872 PetscFunctionReturn(0); 1873 } 1874 1875 #undef __FUNCT__ 1876 #define __FUNCT__ "MatDuplicate_MPIDense" 1877 static PetscErrorCode MatDuplicate_MPIDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 1878 { 1879 Mat mat; 1880 Mat_MPIDense *a,*oldmat = (Mat_MPIDense*)A->data; 1881 PetscErrorCode ierr; 1882 1883 PetscFunctionBegin; 1884 *newmat = 0; 1885 ierr = MatCreate(((PetscObject)A)->comm,&mat);CHKERRQ(ierr); 1886 ierr = MatSetSizes(mat,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1887 ierr = MatSetType(mat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1888 a = (Mat_MPIDense*)mat->data; 1889 ierr = PetscMemcpy(mat->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1890 1891 mat->factortype = A->factortype; 1892 mat->assembled = PETSC_TRUE; 1893 mat->preallocated = PETSC_TRUE; 1894 1895 a->size = oldmat->size; 1896 a->rank = oldmat->rank; 1897 mat->insertmode = NOT_SET_VALUES; 1898 a->nvec = oldmat->nvec; 1899 a->donotstash = oldmat->donotstash; 1900 1901 ierr = PetscLayoutReference(A->rmap,&mat->rmap);CHKERRQ(ierr); 1902 ierr = PetscLayoutReference(A->cmap,&mat->cmap);CHKERRQ(ierr); 1903 1904 ierr = MatSetUpMultiply_MPIDense(mat);CHKERRQ(ierr); 1905 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 1906 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1907 1908 *newmat = mat; 1909 PetscFunctionReturn(0); 1910 } 1911 1912 #undef __FUNCT__ 1913 #define __FUNCT__ "MatLoad_MPIDense_DenseInFile" 1914 PetscErrorCode MatLoad_MPIDense_DenseInFile(MPI_Comm comm,PetscInt fd,PetscInt M,PetscInt N,Mat newmat,PetscInt sizesset) 1915 { 1916 PetscErrorCode ierr; 1917 PetscMPIInt rank,size; 1918 PetscInt *rowners,i,m,nz,j; 1919 PetscScalar *array,*vals,*vals_ptr; 1920 1921 PetscFunctionBegin; 1922 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1923 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1924 1925 /* determine ownership of all rows */ 1926 if (newmat->rmap->n < 0) m = M/size + ((M % size) > rank); 1927 else m = newmat->rmap->n; 1928 ierr = PetscMalloc((size+2)*sizeof(PetscInt),&rowners);CHKERRQ(ierr); 1929 ierr = MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 1930 rowners[0] = 0; 1931 for (i=2; i<=size; i++) { 1932 rowners[i] += rowners[i-1]; 1933 } 1934 1935 if (!sizesset) { 1936 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 1937 } 1938 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1939 ierr = MatGetArray(newmat,&array);CHKERRQ(ierr); 1940 1941 if (!rank) { 1942 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1943 1944 /* read in my part of the matrix numerical values */ 1945 ierr = PetscBinaryRead(fd,vals,m*N,PETSC_SCALAR);CHKERRQ(ierr); 1946 1947 /* insert into matrix-by row (this is why cannot directly read into array */ 1948 vals_ptr = vals; 1949 for (i=0; i<m; i++) { 1950 for (j=0; j<N; j++) { 1951 array[i + j*m] = *vals_ptr++; 1952 } 1953 } 1954 1955 /* read in other processors and ship out */ 1956 for (i=1; i<size; i++) { 1957 nz = (rowners[i+1] - rowners[i])*N; 1958 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 1959 ierr = MPIULong_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1960 } 1961 } else { 1962 /* receive numeric values */ 1963 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1964 1965 /* receive message of values*/ 1966 ierr = MPIULong_Recv(vals,m*N,MPIU_SCALAR,0,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1967 1968 /* insert into matrix-by row (this is why cannot directly read into array */ 1969 vals_ptr = vals; 1970 for (i=0; i<m; i++) { 1971 for (j=0; j<N; j++) { 1972 array[i + j*m] = *vals_ptr++; 1973 } 1974 } 1975 } 1976 ierr = PetscFree(rowners);CHKERRQ(ierr); 1977 ierr = PetscFree(vals);CHKERRQ(ierr); 1978 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1979 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1980 PetscFunctionReturn(0); 1981 } 1982 1983 #undef __FUNCT__ 1984 #define __FUNCT__ "MatLoad_MPIDense" 1985 PetscErrorCode MatLoad_MPIDense(Mat newmat,PetscViewer viewer) 1986 { 1987 PetscScalar *vals,*svals; 1988 MPI_Comm comm = ((PetscObject)viewer)->comm; 1989 MPI_Status status; 1990 PetscMPIInt rank,size,tag = ((PetscObject)viewer)->tag,*rowners,*sndcounts,m,maxnz; 1991 PetscInt header[4],*rowlengths = 0,M,N,*cols; 1992 PetscInt *ourlens,*procsnz = 0,*offlens,jj,*mycols,*smycols; 1993 PetscInt i,nz,j,rstart,rend,sizesset=1,grows,gcols; 1994 int fd; 1995 PetscErrorCode ierr; 1996 1997 PetscFunctionBegin; 1998 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1999 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 2000 if (!rank) { 2001 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2002 ierr = PetscBinaryRead(fd,(char *)header,4,PETSC_INT);CHKERRQ(ierr); 2003 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object"); 2004 } 2005 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) sizesset = 0; 2006 2007 ierr = MPI_Bcast(header+1,3,MPIU_INT,0,comm);CHKERRQ(ierr); 2008 M = header[1]; N = header[2]; nz = header[3]; 2009 2010 /* If global rows/cols are set to PETSC_DECIDE, set it to the sizes given in the file */ 2011 if (sizesset && newmat->rmap->N < 0) newmat->rmap->N = M; 2012 if (sizesset && newmat->cmap->N < 0) newmat->cmap->N = N; 2013 2014 /* If global sizes are set, check if they are consistent with that given in the file */ 2015 if (sizesset) { 2016 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 2017 } 2018 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); 2019 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); 2020 2021 /* 2022 Handle case where matrix is stored on disk as a dense matrix 2023 */ 2024 if (nz == MATRIX_BINARY_FORMAT_DENSE) { 2025 ierr = MatLoad_MPIDense_DenseInFile(comm,fd,M,N,newmat,sizesset);CHKERRQ(ierr); 2026 PetscFunctionReturn(0); 2027 } 2028 2029 /* determine ownership of all rows */ 2030 if (newmat->rmap->n < 0) m = PetscMPIIntCast(M/size + ((M % size) > rank)); 2031 else m = PetscMPIIntCast(newmat->rmap->n); 2032 ierr = PetscMalloc((size+2)*sizeof(PetscMPIInt),&rowners);CHKERRQ(ierr); 2033 ierr = MPI_Allgather(&m,1,MPI_INT,rowners+1,1,MPI_INT,comm);CHKERRQ(ierr); 2034 rowners[0] = 0; 2035 for (i=2; i<=size; i++) { 2036 rowners[i] += rowners[i-1]; 2037 } 2038 rstart = rowners[rank]; 2039 rend = rowners[rank+1]; 2040 2041 /* distribute row lengths to all processors */ 2042 ierr = PetscMalloc2(rend-rstart,PetscInt,&ourlens,rend-rstart,PetscInt,&offlens);CHKERRQ(ierr); 2043 if (!rank) { 2044 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2045 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2046 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&sndcounts);CHKERRQ(ierr); 2047 for (i=0; i<size; i++) sndcounts[i] = rowners[i+1] - rowners[i]; 2048 ierr = MPI_Scatterv(rowlengths,sndcounts,rowners,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2049 ierr = PetscFree(sndcounts);CHKERRQ(ierr); 2050 } else { 2051 ierr = MPI_Scatterv(0,0,0,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2052 } 2053 2054 if (!rank) { 2055 /* calculate the number of nonzeros on each processor */ 2056 ierr = PetscMalloc(size*sizeof(PetscInt),&procsnz);CHKERRQ(ierr); 2057 ierr = PetscMemzero(procsnz,size*sizeof(PetscInt));CHKERRQ(ierr); 2058 for (i=0; i<size; i++) { 2059 for (j=rowners[i]; j< rowners[i+1]; j++) { 2060 procsnz[i] += rowlengths[j]; 2061 } 2062 } 2063 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2064 2065 /* determine max buffer needed and allocate it */ 2066 maxnz = 0; 2067 for (i=0; i<size; i++) { 2068 maxnz = PetscMax(maxnz,procsnz[i]); 2069 } 2070 ierr = PetscMalloc(maxnz*sizeof(PetscInt),&cols);CHKERRQ(ierr); 2071 2072 /* read in my part of the matrix column indices */ 2073 nz = procsnz[0]; 2074 ierr = PetscMalloc(nz*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2075 ierr = PetscBinaryRead(fd,mycols,nz,PETSC_INT);CHKERRQ(ierr); 2076 2077 /* read in every one elses and ship off */ 2078 for (i=1; i<size; i++) { 2079 nz = procsnz[i]; 2080 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 2081 ierr = MPI_Send(cols,nz,MPIU_INT,i,tag,comm);CHKERRQ(ierr); 2082 } 2083 ierr = PetscFree(cols);CHKERRQ(ierr); 2084 } else { 2085 /* determine buffer space needed for message */ 2086 nz = 0; 2087 for (i=0; i<m; i++) { 2088 nz += ourlens[i]; 2089 } 2090 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2091 2092 /* receive message of column indices*/ 2093 ierr = MPI_Recv(mycols,nz,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 2094 ierr = MPI_Get_count(&status,MPIU_INT,&maxnz);CHKERRQ(ierr); 2095 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2096 } 2097 2098 /* loop over local rows, determining number of off diagonal entries */ 2099 ierr = PetscMemzero(offlens,m*sizeof(PetscInt));CHKERRQ(ierr); 2100 jj = 0; 2101 for (i=0; i<m; i++) { 2102 for (j=0; j<ourlens[i]; j++) { 2103 if (mycols[jj] < rstart || mycols[jj] >= rend) offlens[i]++; 2104 jj++; 2105 } 2106 } 2107 2108 /* create our matrix */ 2109 for (i=0; i<m; i++) { 2110 ourlens[i] -= offlens[i]; 2111 } 2112 2113 if (!sizesset) { 2114 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2115 } 2116 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 2117 for (i=0; i<m; i++) { 2118 ourlens[i] += offlens[i]; 2119 } 2120 2121 if (!rank) { 2122 ierr = PetscMalloc(maxnz*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2123 2124 /* read in my part of the matrix numerical values */ 2125 nz = procsnz[0]; 2126 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2127 2128 /* insert into matrix */ 2129 jj = rstart; 2130 smycols = mycols; 2131 svals = vals; 2132 for (i=0; i<m; i++) { 2133 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2134 smycols += ourlens[i]; 2135 svals += ourlens[i]; 2136 jj++; 2137 } 2138 2139 /* read in other processors and ship out */ 2140 for (i=1; i<size; i++) { 2141 nz = procsnz[i]; 2142 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2143 ierr = MPI_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 2144 } 2145 ierr = PetscFree(procsnz);CHKERRQ(ierr); 2146 } else { 2147 /* receive numeric values */ 2148 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2149 2150 /* receive message of values*/ 2151 ierr = MPI_Recv(vals,nz,MPIU_SCALAR,0,((PetscObject)newmat)->tag,comm,&status);CHKERRQ(ierr); 2152 ierr = MPI_Get_count(&status,MPIU_SCALAR,&maxnz);CHKERRQ(ierr); 2153 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2154 2155 /* insert into matrix */ 2156 jj = rstart; 2157 smycols = mycols; 2158 svals = vals; 2159 for (i=0; i<m; i++) { 2160 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2161 smycols += ourlens[i]; 2162 svals += ourlens[i]; 2163 jj++; 2164 } 2165 } 2166 ierr = PetscFree2(ourlens,offlens);CHKERRQ(ierr); 2167 ierr = PetscFree(vals);CHKERRQ(ierr); 2168 ierr = PetscFree(mycols);CHKERRQ(ierr); 2169 ierr = PetscFree(rowners);CHKERRQ(ierr); 2170 2171 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2172 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2173 PetscFunctionReturn(0); 2174 } 2175 2176 #undef __FUNCT__ 2177 #define __FUNCT__ "MatEqual_MPIDense" 2178 PetscErrorCode MatEqual_MPIDense(Mat A,Mat B,PetscBool *flag) 2179 { 2180 Mat_MPIDense *matB = (Mat_MPIDense*)B->data,*matA = (Mat_MPIDense*)A->data; 2181 Mat a,b; 2182 PetscBool flg; 2183 PetscErrorCode ierr; 2184 2185 PetscFunctionBegin; 2186 a = matA->A; 2187 b = matB->A; 2188 ierr = MatEqual(a,b,&flg);CHKERRQ(ierr); 2189 ierr = MPI_Allreduce(&flg,flag,1,MPI_INT,MPI_LAND,((PetscObject)A)->comm);CHKERRQ(ierr); 2190 PetscFunctionReturn(0); 2191 } 2192 2193 #if defined(PETSC_HAVE_PLAPACK) 2194 2195 #undef __FUNCT__ 2196 #define __FUNCT__ "PetscPLAPACKFinalizePackage" 2197 /*@C 2198 PetscPLAPACKFinalizePackage - This function destroys everything in the Petsc interface to PLAPACK. 2199 Level: developer 2200 2201 .keywords: Petsc, destroy, package, PLAPACK 2202 .seealso: PetscFinalize() 2203 @*/ 2204 PetscErrorCode PetscPLAPACKFinalizePackage(void) 2205 { 2206 PetscErrorCode ierr; 2207 2208 PetscFunctionBegin; 2209 ierr = PLA_Finalize();CHKERRQ(ierr); 2210 PetscFunctionReturn(0); 2211 } 2212 2213 #undef __FUNCT__ 2214 #define __FUNCT__ "PetscPLAPACKInitializePackage" 2215 /*@C 2216 PetscPLAPACKInitializePackage - This function initializes everything in the Petsc interface to PLAPACK. It is 2217 called from MatCreate_MPIDense() the first time an MPI dense matrix is called. 2218 2219 Input Parameter: 2220 . comm - the communicator the matrix lives on 2221 2222 Level: developer 2223 2224 Notes: PLAPACK does not have a good fit with MPI communicators; all (parallel) PLAPACK objects have to live in the 2225 same communicator (because there is some global state that is initialized and used for all matrices). In addition if 2226 PLAPACK is initialized (that is the initial matrices created) are on subcommunicators of MPI_COMM_WORLD, these subcommunicators 2227 cannot overlap. 2228 2229 .keywords: Petsc, initialize, package, PLAPACK 2230 .seealso: PetscSysInitializePackage(), PetscInitialize() 2231 @*/ 2232 PetscErrorCode PetscPLAPACKInitializePackage(MPI_Comm comm) 2233 { 2234 PetscMPIInt size; 2235 PetscErrorCode ierr; 2236 2237 PetscFunctionBegin; 2238 if (!PLA_Initialized(PETSC_NULL)) { 2239 2240 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2241 Plapack_nprows = 1; 2242 Plapack_npcols = size; 2243 2244 ierr = PetscOptionsBegin(comm,PETSC_NULL,"PLAPACK Options","Mat");CHKERRQ(ierr); 2245 ierr = PetscOptionsInt("-mat_plapack_nprows","row dimension of 2D processor mesh","None",Plapack_nprows,&Plapack_nprows,PETSC_NULL);CHKERRQ(ierr); 2246 ierr = PetscOptionsInt("-mat_plapack_npcols","column dimension of 2D processor mesh","None",Plapack_npcols,&Plapack_npcols,PETSC_NULL);CHKERRQ(ierr); 2247 #if defined(PETSC_USE_DEBUG) 2248 Plapack_ierror = 3; 2249 #else 2250 Plapack_ierror = 0; 2251 #endif 2252 ierr = PetscOptionsInt("-mat_plapack_ckerror","error checking flag","None",Plapack_ierror,&Plapack_ierror,PETSC_NULL);CHKERRQ(ierr); 2253 if (Plapack_ierror){ 2254 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );CHKERRQ(ierr); 2255 } else { 2256 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );CHKERRQ(ierr); 2257 } 2258 2259 Plapack_nb_alg = 0; 2260 ierr = PetscOptionsInt("-mat_plapack_nb_alg","algorithmic block size","None",Plapack_nb_alg,&Plapack_nb_alg,PETSC_NULL);CHKERRQ(ierr); 2261 if (Plapack_nb_alg) { 2262 ierr = pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,Plapack_nb_alg);CHKERRQ(ierr); 2263 } 2264 PetscOptionsEnd(); 2265 2266 ierr = PLA_Comm_1D_to_2D(comm,Plapack_nprows,Plapack_npcols,&Plapack_comm_2d);CHKERRQ(ierr); 2267 ierr = PLA_Init(Plapack_comm_2d);CHKERRQ(ierr); 2268 ierr = PetscRegisterFinalize(PetscPLAPACKFinalizePackage);CHKERRQ(ierr); 2269 } 2270 PetscFunctionReturn(0); 2271 } 2272 2273 #endif 2274