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_USE_HASH_TABLE: 866 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 867 break; 868 case MAT_IGNORE_OFF_PROC_ENTRIES: 869 a->donotstash = flg; 870 break; 871 case MAT_SYMMETRIC: 872 case MAT_STRUCTURALLY_SYMMETRIC: 873 case MAT_HERMITIAN: 874 case MAT_SYMMETRY_ETERNAL: 875 case MAT_IGNORE_LOWER_TRIANGULAR: 876 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 877 break; 878 default: 879 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 880 } 881 PetscFunctionReturn(0); 882 } 883 884 885 #undef __FUNCT__ 886 #define __FUNCT__ "MatDiagonalScale_MPIDense" 887 PetscErrorCode MatDiagonalScale_MPIDense(Mat A,Vec ll,Vec rr) 888 { 889 Mat_MPIDense *mdn = (Mat_MPIDense*)A->data; 890 Mat_SeqDense *mat = (Mat_SeqDense*)mdn->A->data; 891 PetscScalar *l,*r,x,*v; 892 PetscErrorCode ierr; 893 PetscInt i,j,s2a,s3a,s2,s3,m=mdn->A->rmap->n,n=mdn->A->cmap->n; 894 895 PetscFunctionBegin; 896 ierr = MatGetLocalSize(A,&s2,&s3);CHKERRQ(ierr); 897 if (ll) { 898 ierr = VecGetLocalSize(ll,&s2a);CHKERRQ(ierr); 899 if (s2a != s2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector non-conforming local size, %d != %d.", s2a, s2); 900 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 901 for (i=0; i<m; i++) { 902 x = l[i]; 903 v = mat->v + i; 904 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 905 } 906 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 907 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 908 } 909 if (rr) { 910 ierr = VecGetLocalSize(rr,&s3a);CHKERRQ(ierr); 911 if (s3a != s3) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec non-conforming local size, %d != %d.", s3a, s3); 912 ierr = VecScatterBegin(mdn->Mvctx,rr,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 913 ierr = VecScatterEnd(mdn->Mvctx,rr,mdn->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 914 ierr = VecGetArray(mdn->lvec,&r);CHKERRQ(ierr); 915 for (i=0; i<n; i++) { 916 x = r[i]; 917 v = mat->v + i*m; 918 for (j=0; j<m; j++) { (*v++) *= x;} 919 } 920 ierr = VecRestoreArray(mdn->lvec,&r);CHKERRQ(ierr); 921 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 922 } 923 PetscFunctionReturn(0); 924 } 925 926 #undef __FUNCT__ 927 #define __FUNCT__ "MatNorm_MPIDense" 928 PetscErrorCode MatNorm_MPIDense(Mat A,NormType type,PetscReal *nrm) 929 { 930 Mat_MPIDense *mdn = (Mat_MPIDense*)A->data; 931 Mat_SeqDense *mat = (Mat_SeqDense*)mdn->A->data; 932 PetscErrorCode ierr; 933 PetscInt i,j; 934 PetscReal sum = 0.0; 935 PetscScalar *v = mat->v; 936 937 PetscFunctionBegin; 938 if (mdn->size == 1) { 939 ierr = MatNorm(mdn->A,type,nrm);CHKERRQ(ierr); 940 } else { 941 if (type == NORM_FROBENIUS) { 942 for (i=0; i<mdn->A->cmap->n*mdn->A->rmap->n; i++) { 943 #if defined(PETSC_USE_COMPLEX) 944 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 945 #else 946 sum += (*v)*(*v); v++; 947 #endif 948 } 949 ierr = MPI_Allreduce(&sum,nrm,1,MPIU_REAL,MPIU_SUM,((PetscObject)A)->comm);CHKERRQ(ierr); 950 *nrm = PetscSqrtReal(*nrm); 951 ierr = PetscLogFlops(2.0*mdn->A->cmap->n*mdn->A->rmap->n);CHKERRQ(ierr); 952 } else if (type == NORM_1) { 953 PetscReal *tmp,*tmp2; 954 ierr = PetscMalloc2(A->cmap->N,PetscReal,&tmp,A->cmap->N,PetscReal,&tmp2);CHKERRQ(ierr); 955 ierr = PetscMemzero(tmp,A->cmap->N*sizeof(PetscReal));CHKERRQ(ierr); 956 ierr = PetscMemzero(tmp2,A->cmap->N*sizeof(PetscReal));CHKERRQ(ierr); 957 *nrm = 0.0; 958 v = mat->v; 959 for (j=0; j<mdn->A->cmap->n; j++) { 960 for (i=0; i<mdn->A->rmap->n; i++) { 961 tmp[j] += PetscAbsScalar(*v); v++; 962 } 963 } 964 ierr = MPI_Allreduce(tmp,tmp2,A->cmap->N,MPIU_REAL,MPIU_SUM,((PetscObject)A)->comm);CHKERRQ(ierr); 965 for (j=0; j<A->cmap->N; j++) { 966 if (tmp2[j] > *nrm) *nrm = tmp2[j]; 967 } 968 ierr = PetscFree2(tmp,tmp);CHKERRQ(ierr); 969 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 970 } else if (type == NORM_INFINITY) { /* max row norm */ 971 PetscReal ntemp; 972 ierr = MatNorm(mdn->A,type,&ntemp);CHKERRQ(ierr); 973 ierr = MPI_Allreduce(&ntemp,nrm,1,MPIU_REAL,MPIU_MAX,((PetscObject)A)->comm);CHKERRQ(ierr); 974 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No support for two norm"); 975 } 976 PetscFunctionReturn(0); 977 } 978 979 #undef __FUNCT__ 980 #define __FUNCT__ "MatTranspose_MPIDense" 981 PetscErrorCode MatTranspose_MPIDense(Mat A,MatReuse reuse,Mat *matout) 982 { 983 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 984 Mat_SeqDense *Aloc = (Mat_SeqDense*)a->A->data; 985 Mat B; 986 PetscInt M = A->rmap->N,N = A->cmap->N,m,n,*rwork,rstart = A->rmap->rstart; 987 PetscErrorCode ierr; 988 PetscInt j,i; 989 PetscScalar *v; 990 991 PetscFunctionBegin; 992 if (reuse == MAT_REUSE_MATRIX && A == *matout && M != N) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Supports square matrix only in-place"); 993 if (reuse == MAT_INITIAL_MATRIX || A == *matout) { 994 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 995 ierr = MatSetSizes(B,A->cmap->n,A->rmap->n,N,M);CHKERRQ(ierr); 996 ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr); 997 ierr = MatMPIDenseSetPreallocation(B,PETSC_NULL);CHKERRQ(ierr); 998 } else { 999 B = *matout; 1000 } 1001 1002 m = a->A->rmap->n; n = a->A->cmap->n; v = Aloc->v; 1003 ierr = PetscMalloc(m*sizeof(PetscInt),&rwork);CHKERRQ(ierr); 1004 for (i=0; i<m; i++) rwork[i] = rstart + i; 1005 for (j=0; j<n; j++) { 1006 ierr = MatSetValues(B,1,&j,m,rwork,v,INSERT_VALUES);CHKERRQ(ierr); 1007 v += m; 1008 } 1009 ierr = PetscFree(rwork);CHKERRQ(ierr); 1010 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1012 if (reuse == MAT_INITIAL_MATRIX || *matout != A) { 1013 *matout = B; 1014 } else { 1015 ierr = MatHeaderMerge(A,B);CHKERRQ(ierr); 1016 } 1017 PetscFunctionReturn(0); 1018 } 1019 1020 1021 static PetscErrorCode MatDuplicate_MPIDense(Mat,MatDuplicateOption,Mat *); 1022 extern PetscErrorCode MatScale_MPIDense(Mat,PetscScalar); 1023 1024 #undef __FUNCT__ 1025 #define __FUNCT__ "MatSetUpPreallocation_MPIDense" 1026 PetscErrorCode MatSetUpPreallocation_MPIDense(Mat A) 1027 { 1028 PetscErrorCode ierr; 1029 1030 PetscFunctionBegin; 1031 ierr = MatMPIDenseSetPreallocation(A,0);CHKERRQ(ierr); 1032 PetscFunctionReturn(0); 1033 } 1034 1035 #if defined(PETSC_HAVE_PLAPACK) 1036 1037 #undef __FUNCT__ 1038 #define __FUNCT__ "MatMPIDenseCopyToPlapack" 1039 PetscErrorCode MatMPIDenseCopyToPlapack(Mat A,Mat F) 1040 { 1041 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1042 PetscErrorCode ierr; 1043 PetscInt M=A->cmap->N,m=A->rmap->n,rstart; 1044 PetscScalar *array; 1045 PetscReal one = 1.0; 1046 1047 PetscFunctionBegin; 1048 /* Copy A into F->lu->A */ 1049 ierr = PLA_Obj_set_to_zero(lu->A);CHKERRQ(ierr); 1050 ierr = PLA_API_begin();CHKERRQ(ierr); 1051 ierr = PLA_Obj_API_open(lu->A);CHKERRQ(ierr); 1052 ierr = MatGetOwnershipRange(A,&rstart,PETSC_NULL);CHKERRQ(ierr); 1053 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1054 ierr = PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0);CHKERRQ(ierr); 1055 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1056 ierr = PLA_Obj_API_close(lu->A);CHKERRQ(ierr); 1057 ierr = PLA_API_end();CHKERRQ(ierr); 1058 lu->rstart = rstart; 1059 PetscFunctionReturn(0); 1060 } 1061 1062 #undef __FUNCT__ 1063 #define __FUNCT__ "MatMPIDenseCopyFromPlapack" 1064 PetscErrorCode MatMPIDenseCopyFromPlapack(Mat F,Mat A) 1065 { 1066 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1067 PetscErrorCode ierr; 1068 PetscInt M=A->cmap->N,m=A->rmap->n,rstart; 1069 PetscScalar *array; 1070 PetscReal one = 1.0; 1071 1072 PetscFunctionBegin; 1073 /* Copy F into A->lu->A */ 1074 ierr = MatZeroEntries(A);CHKERRQ(ierr); 1075 ierr = PLA_API_begin();CHKERRQ(ierr); 1076 ierr = PLA_Obj_API_open(lu->A);CHKERRQ(ierr); 1077 ierr = MatGetOwnershipRange(A,&rstart,PETSC_NULL);CHKERRQ(ierr); 1078 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1079 ierr = PLA_API_axpy_global_to_matrix(m,M, &one,lu->A,rstart,0,(void *)array,m);CHKERRQ(ierr); 1080 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1081 ierr = PLA_Obj_API_close(lu->A);CHKERRQ(ierr); 1082 ierr = PLA_API_end();CHKERRQ(ierr); 1083 lu->rstart = rstart; 1084 PetscFunctionReturn(0); 1085 } 1086 1087 #undef __FUNCT__ 1088 #define __FUNCT__ "MatMatMultNumeric_MPIDense_MPIDense" 1089 PetscErrorCode MatMatMultNumeric_MPIDense_MPIDense(Mat A,Mat B,Mat C) 1090 { 1091 PetscErrorCode ierr; 1092 Mat_Plapack *luA = (Mat_Plapack*)A->spptr; 1093 Mat_Plapack *luB = (Mat_Plapack*)B->spptr; 1094 Mat_Plapack *luC = (Mat_Plapack*)C->spptr; 1095 PLA_Obj alpha = NULL,beta = NULL; 1096 1097 PetscFunctionBegin; 1098 ierr = MatMPIDenseCopyToPlapack(A,A);CHKERRQ(ierr); 1099 ierr = MatMPIDenseCopyToPlapack(B,B);CHKERRQ(ierr); 1100 1101 /* 1102 ierr = PLA_Global_show("A = ",luA->A,"%g ","");CHKERRQ(ierr); 1103 ierr = PLA_Global_show("B = ",luB->A,"%g ","");CHKERRQ(ierr); 1104 */ 1105 1106 /* do the multiply in PLA */ 1107 ierr = PLA_Create_constants_conf_to(luA->A,NULL,NULL,&alpha);CHKERRQ(ierr); 1108 ierr = PLA_Create_constants_conf_to(luC->A,NULL,&beta,NULL);CHKERRQ(ierr); 1109 CHKMEMQ; 1110 1111 ierr = PLA_Gemm(PLA_NO_TRANSPOSE,PLA_NO_TRANSPOSE,alpha,luA->A,luB->A,beta,luC->A); /* CHKERRQ(ierr); */ 1112 CHKMEMQ; 1113 ierr = PLA_Obj_free(&alpha);CHKERRQ(ierr); 1114 ierr = PLA_Obj_free(&beta);CHKERRQ(ierr); 1115 1116 /* 1117 ierr = PLA_Global_show("C = ",luC->A,"%g ","");CHKERRQ(ierr); 1118 */ 1119 ierr = MatMPIDenseCopyFromPlapack(C,C);CHKERRQ(ierr); 1120 PetscFunctionReturn(0); 1121 } 1122 1123 #undef __FUNCT__ 1124 #define __FUNCT__ "MatMatMultSymbolic_MPIDense_MPIDense" 1125 PetscErrorCode MatMatMultSymbolic_MPIDense_MPIDense(Mat A,Mat B,PetscReal fill,Mat *C) 1126 { 1127 PetscErrorCode ierr; 1128 PetscInt m=A->rmap->n,n=B->cmap->n; 1129 Mat Cmat; 1130 1131 PetscFunctionBegin; 1132 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); 1133 SETERRQ(((PetscObject)A)->comm,PETSC_ERR_LIB,"Due to apparent bugs in PLAPACK,this is not currently supported"); 1134 ierr = MatCreate(((PetscObject)B)->comm,&Cmat);CHKERRQ(ierr); 1135 ierr = MatSetSizes(Cmat,m,n,A->rmap->N,B->cmap->N);CHKERRQ(ierr); 1136 ierr = MatSetType(Cmat,MATMPIDENSE);CHKERRQ(ierr); 1137 ierr = MatAssemblyBegin(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1138 ierr = MatAssemblyEnd(Cmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1139 1140 *C = Cmat; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "MatMatMult_MPIDense_MPIDense" 1146 PetscErrorCode MatMatMult_MPIDense_MPIDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1147 { 1148 PetscErrorCode ierr; 1149 1150 PetscFunctionBegin; 1151 if (scall == MAT_INITIAL_MATRIX){ 1152 ierr = MatMatMultSymbolic_MPIDense_MPIDense(A,B,fill,C);CHKERRQ(ierr); 1153 } 1154 ierr = MatMatMultNumeric_MPIDense_MPIDense(A,B,*C);CHKERRQ(ierr); 1155 PetscFunctionReturn(0); 1156 } 1157 1158 #undef __FUNCT__ 1159 #define __FUNCT__ "MatSolve_MPIDense" 1160 PetscErrorCode MatSolve_MPIDense(Mat A,Vec b,Vec x) 1161 { 1162 MPI_Comm comm = ((PetscObject)A)->comm; 1163 Mat_Plapack *lu = (Mat_Plapack*)A->spptr; 1164 PetscErrorCode ierr; 1165 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,i,j,*idx_pla,*idx_petsc,loc_m,loc_stride; 1166 PetscScalar *array; 1167 PetscReal one = 1.0; 1168 PetscMPIInt size,rank,r_rank,r_nproc,c_rank,c_nproc;; 1169 PLA_Obj v_pla = NULL; 1170 PetscScalar *loc_buf; 1171 Vec loc_x; 1172 1173 PetscFunctionBegin; 1174 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1175 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1176 1177 /* Create PLAPACK vector objects, then copy b into PLAPACK b */ 1178 PLA_Mvector_create(lu->datatype,M,1,lu->templ,PLA_ALIGN_FIRST,&v_pla); 1179 PLA_Obj_set_to_zero(v_pla); 1180 1181 /* Copy b into rhs_pla */ 1182 PLA_API_begin(); 1183 PLA_Obj_API_open(v_pla); 1184 ierr = VecGetArray(b,&array);CHKERRQ(ierr); 1185 PLA_API_axpy_vector_to_global(m,&one,(void *)array,1,v_pla,lu->rstart); 1186 ierr = VecRestoreArray(b,&array);CHKERRQ(ierr); 1187 PLA_Obj_API_close(v_pla); 1188 PLA_API_end(); 1189 1190 if (A->factortype == MAT_FACTOR_LU){ 1191 /* Apply the permutations to the right hand sides */ 1192 PLA_Apply_pivots_to_rows (v_pla,lu->pivots); 1193 1194 /* Solve L y = b, overwriting b with y */ 1195 PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_UNIT_DIAG,lu->A,v_pla ); 1196 1197 /* Solve U x = y (=b), overwriting b with x */ 1198 PLA_Trsv( PLA_UPPER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla ); 1199 } else { /* MAT_FACTOR_CHOLESKY */ 1200 PLA_Trsv( PLA_LOWER_TRIANGULAR,PLA_NO_TRANSPOSE,PLA_NONUNIT_DIAG,lu->A,v_pla); 1201 PLA_Trsv( PLA_LOWER_TRIANGULAR,(lu->datatype == MPI_DOUBLE ? PLA_TRANSPOSE : PLA_CONJUGATE_TRANSPOSE), 1202 PLA_NONUNIT_DIAG,lu->A,v_pla); 1203 } 1204 1205 /* Copy PLAPACK x into Petsc vector x */ 1206 PLA_Obj_local_length(v_pla, &loc_m); 1207 PLA_Obj_local_buffer(v_pla, (void**)&loc_buf); 1208 PLA_Obj_local_stride(v_pla, &loc_stride); 1209 /* 1210 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); 1211 */ 1212 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,loc_m*loc_stride,loc_buf,&loc_x);CHKERRQ(ierr); 1213 if (!lu->pla_solved){ 1214 1215 PLA_Temp_comm_row_info(lu->templ,&Plapack_comm_2d,&r_rank,&r_nproc); 1216 PLA_Temp_comm_col_info(lu->templ,&Plapack_comm_2d,&c_rank,&c_nproc); 1217 1218 /* Create IS and cts for VecScatterring */ 1219 PLA_Obj_local_length(v_pla, &loc_m); 1220 PLA_Obj_local_stride(v_pla, &loc_stride); 1221 ierr = PetscMalloc2(loc_m,PetscInt,&idx_pla,loc_m,PetscInt,&idx_petsc);CHKERRQ(ierr); 1222 1223 rstart = (r_rank*c_nproc+c_rank)*lu->nb; 1224 for (i=0; i<loc_m; i+=lu->nb){ 1225 j = 0; 1226 while (j < lu->nb && i+j < loc_m){ 1227 idx_petsc[i+j] = rstart + j; j++; 1228 } 1229 rstart += size*lu->nb; 1230 } 1231 1232 for (i=0; i<loc_m; i++) idx_pla[i] = i*loc_stride; 1233 1234 ierr = ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_pla,PETSC_COPY_VALUES,&lu->is_pla);CHKERRQ(ierr); 1235 ierr = ISCreateGeneral(PETSC_COMM_SELF,loc_m,idx_petsc,PETSC_COPY_VALUES,&lu->is_petsc);CHKERRQ(ierr); 1236 ierr = PetscFree2(idx_pla,idx_petsc);CHKERRQ(ierr); 1237 ierr = VecScatterCreate(loc_x,lu->is_pla,x,lu->is_petsc,&lu->ctx);CHKERRQ(ierr); 1238 } 1239 ierr = VecScatterBegin(lu->ctx,loc_x,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1240 ierr = VecScatterEnd(lu->ctx,loc_x,x,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1241 1242 /* Free data */ 1243 ierr = VecDestroy(&loc_x);CHKERRQ(ierr); 1244 PLA_Obj_free(&v_pla); 1245 1246 lu->pla_solved = PETSC_TRUE; 1247 PetscFunctionReturn(0); 1248 } 1249 1250 #undef __FUNCT__ 1251 #define __FUNCT__ "MatLUFactorNumeric_MPIDense" 1252 PetscErrorCode MatLUFactorNumeric_MPIDense(Mat F,Mat A,const MatFactorInfo *info) 1253 { 1254 Mat_Plapack *lu = (Mat_Plapack*)(F)->spptr; 1255 PetscErrorCode ierr; 1256 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,rend; 1257 PetscInt info_pla=0; 1258 PetscScalar *array,one = 1.0; 1259 1260 PetscFunctionBegin; 1261 if (lu->mstruct == SAME_NONZERO_PATTERN){ 1262 PLA_Obj_free(&lu->A); 1263 PLA_Obj_free (&lu->pivots); 1264 } 1265 /* Create PLAPACK matrix object */ 1266 lu->A = NULL; lu->pivots = NULL; 1267 PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A); 1268 PLA_Obj_set_to_zero(lu->A); 1269 PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots); 1270 1271 /* Copy A into lu->A */ 1272 PLA_API_begin(); 1273 PLA_Obj_API_open(lu->A); 1274 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1275 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1276 PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0); 1277 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1278 PLA_Obj_API_close(lu->A); 1279 PLA_API_end(); 1280 1281 /* Factor P A -> L U overwriting lower triangular portion of A with L, upper, U */ 1282 info_pla = PLA_LU(lu->A,lu->pivots); 1283 if (info_pla != 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot encountered at row %d from PLA_LU()",info_pla); 1284 1285 lu->rstart = rstart; 1286 lu->mstruct = SAME_NONZERO_PATTERN; 1287 F->ops->solve = MatSolve_MPIDense; 1288 F->assembled = PETSC_TRUE; /* required by -ksp_view */ 1289 PetscFunctionReturn(0); 1290 } 1291 1292 #undef __FUNCT__ 1293 #define __FUNCT__ "MatCholeskyFactorNumeric_MPIDense" 1294 PetscErrorCode MatCholeskyFactorNumeric_MPIDense(Mat F,Mat A,const MatFactorInfo *info) 1295 { 1296 Mat_Plapack *lu = (Mat_Plapack*)F->spptr; 1297 PetscErrorCode ierr; 1298 PetscInt M=A->rmap->N,m=A->rmap->n,rstart,rend; 1299 PetscInt info_pla=0; 1300 PetscScalar *array,one = 1.0; 1301 1302 PetscFunctionBegin; 1303 if (lu->mstruct == SAME_NONZERO_PATTERN){ 1304 PLA_Obj_free(&lu->A); 1305 } 1306 /* Create PLAPACK matrix object */ 1307 lu->A = NULL; 1308 lu->pivots = NULL; 1309 PLA_Matrix_create(lu->datatype,M,M,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A); 1310 1311 /* Copy A into lu->A */ 1312 PLA_API_begin(); 1313 PLA_Obj_API_open(lu->A); 1314 ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 1315 ierr = MatGetArray(A,&array);CHKERRQ(ierr); 1316 PLA_API_axpy_matrix_to_global(m,M, &one,(void *)array,m,lu->A,rstart,0); 1317 ierr = MatRestoreArray(A,&array);CHKERRQ(ierr); 1318 PLA_Obj_API_close(lu->A); 1319 PLA_API_end(); 1320 1321 /* Factor P A -> Chol */ 1322 info_pla = PLA_Chol(PLA_LOWER_TRIANGULAR,lu->A); 1323 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); 1324 1325 lu->rstart = rstart; 1326 lu->mstruct = SAME_NONZERO_PATTERN; 1327 F->ops->solve = MatSolve_MPIDense; 1328 F->assembled = PETSC_TRUE; /* required by -ksp_view */ 1329 PetscFunctionReturn(0); 1330 } 1331 1332 /* Note the Petsc perm permutation is ignored */ 1333 #undef __FUNCT__ 1334 #define __FUNCT__ "MatCholeskyFactorSymbolic_MPIDense" 1335 PetscErrorCode MatCholeskyFactorSymbolic_MPIDense(Mat F,Mat A,IS perm,const MatFactorInfo *info) 1336 { 1337 PetscErrorCode ierr; 1338 PetscBool issymmetric,set; 1339 1340 PetscFunctionBegin; 1341 ierr = MatIsSymmetricKnown(A,&set,&issymmetric);CHKERRQ(ierr); 1342 if (!set || !issymmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be set as MAT_SYMMETRIC for CholeskyFactor()"); 1343 F->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_MPIDense; 1344 PetscFunctionReturn(0); 1345 } 1346 1347 /* Note the Petsc r and c permutations are ignored */ 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "MatLUFactorSymbolic_MPIDense" 1350 PetscErrorCode MatLUFactorSymbolic_MPIDense(Mat F,Mat A,IS r,IS c,const MatFactorInfo *info) 1351 { 1352 PetscErrorCode ierr; 1353 PetscInt M = A->rmap->N; 1354 Mat_Plapack *lu; 1355 1356 PetscFunctionBegin; 1357 lu = (Mat_Plapack*)F->spptr; 1358 ierr = PLA_Mvector_create(MPI_INT,M,1,lu->templ,PLA_ALIGN_FIRST,&lu->pivots);CHKERRQ(ierr); 1359 F->ops->lufactornumeric = MatLUFactorNumeric_MPIDense; 1360 PetscFunctionReturn(0); 1361 } 1362 1363 EXTERN_C_BEGIN 1364 #undef __FUNCT__ 1365 #define __FUNCT__ "MatFactorGetSolverPackage_mpidense_plapack" 1366 PetscErrorCode MatFactorGetSolverPackage_mpidense_plapack(Mat A,const MatSolverPackage *type) 1367 { 1368 PetscFunctionBegin; 1369 *type = MATSOLVERPLAPACK; 1370 PetscFunctionReturn(0); 1371 } 1372 EXTERN_C_END 1373 1374 EXTERN_C_BEGIN 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "MatGetFactor_mpidense_plapack" 1377 PetscErrorCode MatGetFactor_mpidense_plapack(Mat A,MatFactorType ftype,Mat *F) 1378 { 1379 PetscErrorCode ierr; 1380 Mat_Plapack *lu; 1381 PetscMPIInt size; 1382 PetscInt M=A->rmap->N; 1383 1384 PetscFunctionBegin; 1385 /* Create the factorization matrix */ 1386 ierr = MatCreate(((PetscObject)A)->comm,F);CHKERRQ(ierr); 1387 ierr = MatSetSizes(*F,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1388 ierr = MatSetType(*F,((PetscObject)A)->type_name);CHKERRQ(ierr); 1389 ierr = PetscNewLog(*F,Mat_Plapack,&lu);CHKERRQ(ierr); 1390 (*F)->spptr = (void*)lu; 1391 1392 /* Set default Plapack parameters */ 1393 ierr = MPI_Comm_size(((PetscObject)A)->comm,&size);CHKERRQ(ierr); 1394 lu->nb = M/size; 1395 if (M - lu->nb*size) lu->nb++; /* without cyclic distribution */ 1396 1397 /* Set runtime options */ 1398 ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"PLAPACK Options","Mat");CHKERRQ(ierr); 1399 ierr = PetscOptionsInt("-mat_plapack_nb","block size of template vector","None",lu->nb,&lu->nb,PETSC_NULL);CHKERRQ(ierr); 1400 PetscOptionsEnd(); 1401 1402 /* Create object distribution template */ 1403 lu->templ = NULL; 1404 ierr = PLA_Temp_create(lu->nb, 0, &lu->templ);CHKERRQ(ierr); 1405 1406 /* Set the datatype */ 1407 #if defined(PETSC_USE_COMPLEX) 1408 lu->datatype = MPI_DOUBLE_COMPLEX; 1409 #else 1410 lu->datatype = MPI_DOUBLE; 1411 #endif 1412 1413 ierr = PLA_Matrix_create(lu->datatype,M,A->cmap->N,lu->templ,PLA_ALIGN_FIRST,PLA_ALIGN_FIRST,&lu->A);CHKERRQ(ierr); 1414 1415 1416 lu->pla_solved = PETSC_FALSE; /* MatSolve_Plapack() is called yet */ 1417 lu->mstruct = DIFFERENT_NONZERO_PATTERN; 1418 1419 if (ftype == MAT_FACTOR_LU) { 1420 (*F)->ops->lufactorsymbolic = MatLUFactorSymbolic_MPIDense; 1421 } else if (ftype == MAT_FACTOR_CHOLESKY) { 1422 (*F)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_MPIDense; 1423 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"No incomplete factorizations for dense matrices"); 1424 (*F)->factortype = ftype; 1425 ierr = PetscObjectComposeFunctionDynamic((PetscObject)(*F),"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_mpidense_plapack",MatFactorGetSolverPackage_mpidense_plapack);CHKERRQ(ierr); 1426 PetscFunctionReturn(0); 1427 } 1428 EXTERN_C_END 1429 #endif 1430 1431 #undef __FUNCT__ 1432 #define __FUNCT__ "MatGetFactor_mpidense_petsc" 1433 PetscErrorCode MatGetFactor_mpidense_petsc(Mat A,MatFactorType ftype,Mat *F) 1434 { 1435 #if defined(PETSC_HAVE_PLAPACK) 1436 PetscErrorCode ierr; 1437 #endif 1438 1439 PetscFunctionBegin; 1440 #if defined(PETSC_HAVE_PLAPACK) 1441 ierr = MatGetFactor_mpidense_plapack(A,ftype,F);CHKERRQ(ierr); 1442 #else 1443 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix format %s uses PLAPACK direct solver. Install PLAPACK",((PetscObject)A)->type_name); 1444 #endif 1445 PetscFunctionReturn(0); 1446 } 1447 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "MatAXPY_MPIDense" 1450 PetscErrorCode MatAXPY_MPIDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 1451 { 1452 PetscErrorCode ierr; 1453 Mat_MPIDense *A = (Mat_MPIDense*)Y->data, *B = (Mat_MPIDense*)X->data; 1454 1455 PetscFunctionBegin; 1456 ierr = MatAXPY(A->A,alpha,B->A,str);CHKERRQ(ierr); 1457 PetscFunctionReturn(0); 1458 } 1459 1460 #undef __FUNCT__ 1461 #define __FUNCT__ "MatConjugate_MPIDense" 1462 PetscErrorCode MatConjugate_MPIDense(Mat mat) 1463 { 1464 Mat_MPIDense *a = (Mat_MPIDense *)mat->data; 1465 PetscErrorCode ierr; 1466 1467 PetscFunctionBegin; 1468 ierr = MatConjugate(a->A);CHKERRQ(ierr); 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "MatRealPart_MPIDense" 1474 PetscErrorCode MatRealPart_MPIDense(Mat A) 1475 { 1476 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1477 PetscErrorCode ierr; 1478 1479 PetscFunctionBegin; 1480 ierr = MatRealPart(a->A);CHKERRQ(ierr); 1481 PetscFunctionReturn(0); 1482 } 1483 1484 #undef __FUNCT__ 1485 #define __FUNCT__ "MatImaginaryPart_MPIDense" 1486 PetscErrorCode MatImaginaryPart_MPIDense(Mat A) 1487 { 1488 Mat_MPIDense *a = (Mat_MPIDense*)A->data; 1489 PetscErrorCode ierr; 1490 1491 PetscFunctionBegin; 1492 ierr = MatImaginaryPart(a->A);CHKERRQ(ierr); 1493 PetscFunctionReturn(0); 1494 } 1495 1496 extern PetscErrorCode MatGetColumnNorms_SeqDense(Mat,NormType,PetscReal*); 1497 #undef __FUNCT__ 1498 #define __FUNCT__ "MatGetColumnNorms_MPIDense" 1499 PetscErrorCode MatGetColumnNorms_MPIDense(Mat A,NormType type,PetscReal *norms) 1500 { 1501 PetscErrorCode ierr; 1502 PetscInt i,n; 1503 Mat_MPIDense *a = (Mat_MPIDense*) A->data; 1504 PetscReal *work; 1505 1506 PetscFunctionBegin; 1507 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 1508 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 1509 ierr = MatGetColumnNorms_SeqDense(a->A,type,work);CHKERRQ(ierr); 1510 if (type == NORM_2) { 1511 for (i=0; i<n; i++) work[i] *= work[i]; 1512 } 1513 if (type == NORM_INFINITY) { 1514 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_MAX,A->hdr.comm);CHKERRQ(ierr); 1515 } else { 1516 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPIU_SUM,A->hdr.comm);CHKERRQ(ierr); 1517 } 1518 ierr = PetscFree(work);CHKERRQ(ierr); 1519 if (type == NORM_2) { 1520 for (i=0; i<n; i++) norms[i] = PetscSqrtReal(norms[i]); 1521 } 1522 PetscFunctionReturn(0); 1523 } 1524 1525 /* -------------------------------------------------------------------*/ 1526 static struct _MatOps MatOps_Values = {MatSetValues_MPIDense, 1527 MatGetRow_MPIDense, 1528 MatRestoreRow_MPIDense, 1529 MatMult_MPIDense, 1530 /* 4*/ MatMultAdd_MPIDense, 1531 MatMultTranspose_MPIDense, 1532 MatMultTransposeAdd_MPIDense, 1533 0, 1534 0, 1535 0, 1536 /*10*/ 0, 1537 0, 1538 0, 1539 0, 1540 MatTranspose_MPIDense, 1541 /*15*/ MatGetInfo_MPIDense, 1542 MatEqual_MPIDense, 1543 MatGetDiagonal_MPIDense, 1544 MatDiagonalScale_MPIDense, 1545 MatNorm_MPIDense, 1546 /*20*/ MatAssemblyBegin_MPIDense, 1547 MatAssemblyEnd_MPIDense, 1548 MatSetOption_MPIDense, 1549 MatZeroEntries_MPIDense, 1550 /*24*/ MatZeroRows_MPIDense, 1551 0, 1552 0, 1553 0, 1554 0, 1555 /*29*/ MatSetUpPreallocation_MPIDense, 1556 0, 1557 0, 1558 MatGetArray_MPIDense, 1559 MatRestoreArray_MPIDense, 1560 /*34*/ MatDuplicate_MPIDense, 1561 0, 1562 0, 1563 0, 1564 0, 1565 /*39*/ MatAXPY_MPIDense, 1566 MatGetSubMatrices_MPIDense, 1567 0, 1568 MatGetValues_MPIDense, 1569 0, 1570 /*44*/ 0, 1571 MatScale_MPIDense, 1572 0, 1573 0, 1574 0, 1575 /*49*/ 0, 1576 0, 1577 0, 1578 0, 1579 0, 1580 /*54*/ 0, 1581 0, 1582 0, 1583 0, 1584 0, 1585 /*59*/ MatGetSubMatrix_MPIDense, 1586 MatDestroy_MPIDense, 1587 MatView_MPIDense, 1588 0, 1589 0, 1590 /*64*/ 0, 1591 0, 1592 0, 1593 0, 1594 0, 1595 /*69*/ 0, 1596 0, 1597 0, 1598 0, 1599 0, 1600 /*74*/ 0, 1601 0, 1602 0, 1603 0, 1604 0, 1605 /*79*/ 0, 1606 0, 1607 0, 1608 0, 1609 /*83*/ MatLoad_MPIDense, 1610 0, 1611 0, 1612 0, 1613 0, 1614 0, 1615 /*89*/ 1616 #if defined(PETSC_HAVE_PLAPACK) 1617 MatMatMult_MPIDense_MPIDense, 1618 MatMatMultSymbolic_MPIDense_MPIDense, 1619 MatMatMultNumeric_MPIDense_MPIDense, 1620 #else 1621 0, 1622 0, 1623 0, 1624 #endif 1625 0, 1626 0, 1627 /*94*/ 0, 1628 0, 1629 0, 1630 0, 1631 0, 1632 /*99*/ 0, 1633 0, 1634 0, 1635 MatConjugate_MPIDense, 1636 0, 1637 /*104*/0, 1638 MatRealPart_MPIDense, 1639 MatImaginaryPart_MPIDense, 1640 0, 1641 0, 1642 /*109*/0, 1643 0, 1644 0, 1645 0, 1646 0, 1647 /*114*/0, 1648 0, 1649 0, 1650 0, 1651 0, 1652 /*119*/0, 1653 0, 1654 0, 1655 0, 1656 0, 1657 /*124*/0, 1658 MatGetColumnNorms_MPIDense 1659 }; 1660 1661 EXTERN_C_BEGIN 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "MatMPIDenseSetPreallocation_MPIDense" 1664 PetscErrorCode MatMPIDenseSetPreallocation_MPIDense(Mat mat,PetscScalar *data) 1665 { 1666 Mat_MPIDense *a; 1667 PetscErrorCode ierr; 1668 1669 PetscFunctionBegin; 1670 mat->preallocated = PETSC_TRUE; 1671 /* Note: For now, when data is specified above, this assumes the user correctly 1672 allocates the local dense storage space. We should add error checking. */ 1673 1674 a = (Mat_MPIDense*)mat->data; 1675 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 1676 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 1677 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 1678 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 1679 a->nvec = mat->cmap->n; 1680 1681 ierr = MatCreate(PETSC_COMM_SELF,&a->A);CHKERRQ(ierr); 1682 ierr = MatSetSizes(a->A,mat->rmap->n,mat->cmap->N,mat->rmap->n,mat->cmap->N);CHKERRQ(ierr); 1683 ierr = MatSetType(a->A,MATSEQDENSE);CHKERRQ(ierr); 1684 ierr = MatSeqDenseSetPreallocation(a->A,data);CHKERRQ(ierr); 1685 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1686 PetscFunctionReturn(0); 1687 } 1688 EXTERN_C_END 1689 1690 /*MC 1691 MATSOLVERPLAPACK = "mpidense" - Parallel LU and Cholesky factorization for MATMPIDENSE matrices 1692 1693 run ./configure with the option --download-plapack 1694 1695 1696 Options Database Keys: 1697 . -mat_plapack_nprows <n> - number of rows in processor partition 1698 . -mat_plapack_npcols <n> - number of columns in processor partition 1699 . -mat_plapack_nb <n> - block size of template vector 1700 . -mat_plapack_nb_alg <n> - algorithmic block size 1701 - -mat_plapack_ckerror <n> - error checking flag 1702 1703 Level: intermediate 1704 1705 .seealso: MatCreateMPIDense(), MATDENSE, MATSEQDENSE, PCFactorSetSolverPackage(), MatSolverPackage 1706 1707 M*/ 1708 1709 EXTERN_C_BEGIN 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "MatCreate_MPIDense" 1712 PetscErrorCode MatCreate_MPIDense(Mat mat) 1713 { 1714 Mat_MPIDense *a; 1715 PetscErrorCode ierr; 1716 1717 PetscFunctionBegin; 1718 ierr = PetscNewLog(mat,Mat_MPIDense,&a);CHKERRQ(ierr); 1719 mat->data = (void*)a; 1720 ierr = PetscMemcpy(mat->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1721 1722 mat->insertmode = NOT_SET_VALUES; 1723 ierr = MPI_Comm_rank(((PetscObject)mat)->comm,&a->rank);CHKERRQ(ierr); 1724 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&a->size);CHKERRQ(ierr); 1725 1726 /* build cache for off array entries formed */ 1727 a->donotstash = PETSC_FALSE; 1728 ierr = MatStashCreate_Private(((PetscObject)mat)->comm,1,&mat->stash);CHKERRQ(ierr); 1729 1730 /* stuff used for matrix vector multiply */ 1731 a->lvec = 0; 1732 a->Mvctx = 0; 1733 a->roworiented = PETSC_TRUE; 1734 1735 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetDiagonalBlock_C", 1736 "MatGetDiagonalBlock_MPIDense", 1737 MatGetDiagonalBlock_MPIDense);CHKERRQ(ierr); 1738 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMPIDenseSetPreallocation_C", 1739 "MatMPIDenseSetPreallocation_MPIDense", 1740 MatMPIDenseSetPreallocation_MPIDense);CHKERRQ(ierr); 1741 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_mpiaij_mpidense_C", 1742 "MatMatMult_MPIAIJ_MPIDense", 1743 MatMatMult_MPIAIJ_MPIDense);CHKERRQ(ierr); 1744 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_mpiaij_mpidense_C", 1745 "MatMatMultSymbolic_MPIAIJ_MPIDense", 1746 MatMatMultSymbolic_MPIAIJ_MPIDense);CHKERRQ(ierr); 1747 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_mpiaij_mpidense_C", 1748 "MatMatMultNumeric_MPIAIJ_MPIDense", 1749 MatMatMultNumeric_MPIAIJ_MPIDense);CHKERRQ(ierr); 1750 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_petsc_C", 1751 "MatGetFactor_mpidense_petsc", 1752 MatGetFactor_mpidense_petsc);CHKERRQ(ierr); 1753 #if defined(PETSC_HAVE_PLAPACK) 1754 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatGetFactor_plapack_C", 1755 "MatGetFactor_mpidense_plapack", 1756 MatGetFactor_mpidense_plapack);CHKERRQ(ierr); 1757 ierr = PetscPLAPACKInitializePackage(((PetscObject)mat)->comm);CHKERRQ(ierr); 1758 #endif 1759 ierr = PetscObjectChangeTypeName((PetscObject)mat,MATMPIDENSE);CHKERRQ(ierr); 1760 1761 PetscFunctionReturn(0); 1762 } 1763 EXTERN_C_END 1764 1765 /*MC 1766 MATDENSE - MATDENSE = "dense" - A matrix type to be used for dense matrices. 1767 1768 This matrix type is identical to MATSEQDENSE when constructed with a single process communicator, 1769 and MATMPIDENSE otherwise. 1770 1771 Options Database Keys: 1772 . -mat_type dense - sets the matrix type to "dense" during a call to MatSetFromOptions() 1773 1774 Level: beginner 1775 1776 1777 .seealso: MatCreateMPIDense,MATSEQDENSE,MATMPIDENSE 1778 M*/ 1779 1780 #undef __FUNCT__ 1781 #define __FUNCT__ "MatMPIDenseSetPreallocation" 1782 /*@C 1783 MatMPIDenseSetPreallocation - Sets the array used to store the matrix entries 1784 1785 Not collective 1786 1787 Input Parameters: 1788 . A - the matrix 1789 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1790 to control all matrix memory allocation. 1791 1792 Notes: 1793 The dense format is fully compatible with standard Fortran 77 1794 storage by columns. 1795 1796 The data input variable is intended primarily for Fortran programmers 1797 who wish to allocate their own matrix memory space. Most users should 1798 set data=PETSC_NULL. 1799 1800 Level: intermediate 1801 1802 .keywords: matrix,dense, parallel 1803 1804 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1805 @*/ 1806 PetscErrorCode MatMPIDenseSetPreallocation(Mat mat,PetscScalar *data) 1807 { 1808 PetscErrorCode ierr; 1809 1810 PetscFunctionBegin; 1811 ierr = PetscTryMethod(mat,"MatMPIDenseSetPreallocation_C",(Mat,PetscScalar *),(mat,data));CHKERRQ(ierr); 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "MatCreateMPIDense" 1817 /*@C 1818 MatCreateMPIDense - Creates a sparse parallel matrix in dense format. 1819 1820 Collective on MPI_Comm 1821 1822 Input Parameters: 1823 + comm - MPI communicator 1824 . m - number of local rows (or PETSC_DECIDE to have calculated if M is given) 1825 . n - number of local columns (or PETSC_DECIDE to have calculated if N is given) 1826 . M - number of global rows (or PETSC_DECIDE to have calculated if m is given) 1827 . N - number of global columns (or PETSC_DECIDE to have calculated if n is given) 1828 - data - optional location of matrix data. Set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users) for PETSc 1829 to control all matrix memory allocation. 1830 1831 Output Parameter: 1832 . A - the matrix 1833 1834 Notes: 1835 The dense format is fully compatible with standard Fortran 77 1836 storage by columns. 1837 1838 The data input variable is intended primarily for Fortran programmers 1839 who wish to allocate their own matrix memory space. Most users should 1840 set data=PETSC_NULL (PETSC_NULL_SCALAR for Fortran users). 1841 1842 The user MUST specify either the local or global matrix dimensions 1843 (possibly both). 1844 1845 Level: intermediate 1846 1847 .keywords: matrix,dense, parallel 1848 1849 .seealso: MatCreate(), MatCreateSeqDense(), MatSetValues() 1850 @*/ 1851 PetscErrorCode MatCreateMPIDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscScalar *data,Mat *A) 1852 { 1853 PetscErrorCode ierr; 1854 PetscMPIInt size; 1855 1856 PetscFunctionBegin; 1857 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1858 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 1859 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1860 if (size > 1) { 1861 ierr = MatSetType(*A,MATMPIDENSE);CHKERRQ(ierr); 1862 ierr = MatMPIDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1863 } else { 1864 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1865 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1866 } 1867 PetscFunctionReturn(0); 1868 } 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatDuplicate_MPIDense" 1872 static PetscErrorCode MatDuplicate_MPIDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 1873 { 1874 Mat mat; 1875 Mat_MPIDense *a,*oldmat = (Mat_MPIDense*)A->data; 1876 PetscErrorCode ierr; 1877 1878 PetscFunctionBegin; 1879 *newmat = 0; 1880 ierr = MatCreate(((PetscObject)A)->comm,&mat);CHKERRQ(ierr); 1881 ierr = MatSetSizes(mat,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);CHKERRQ(ierr); 1882 ierr = MatSetType(mat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1883 a = (Mat_MPIDense*)mat->data; 1884 ierr = PetscMemcpy(mat->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1885 1886 mat->factortype = A->factortype; 1887 mat->assembled = PETSC_TRUE; 1888 mat->preallocated = PETSC_TRUE; 1889 1890 a->size = oldmat->size; 1891 a->rank = oldmat->rank; 1892 mat->insertmode = NOT_SET_VALUES; 1893 a->nvec = oldmat->nvec; 1894 a->donotstash = oldmat->donotstash; 1895 1896 ierr = PetscLayoutReference(A->rmap,&mat->rmap);CHKERRQ(ierr); 1897 ierr = PetscLayoutReference(A->cmap,&mat->cmap);CHKERRQ(ierr); 1898 1899 ierr = MatSetUpMultiply_MPIDense(mat);CHKERRQ(ierr); 1900 ierr = MatDuplicate(oldmat->A,cpvalues,&a->A);CHKERRQ(ierr); 1901 ierr = PetscLogObjectParent(mat,a->A);CHKERRQ(ierr); 1902 1903 *newmat = mat; 1904 PetscFunctionReturn(0); 1905 } 1906 1907 #undef __FUNCT__ 1908 #define __FUNCT__ "MatLoad_MPIDense_DenseInFile" 1909 PetscErrorCode MatLoad_MPIDense_DenseInFile(MPI_Comm comm,PetscInt fd,PetscInt M,PetscInt N,Mat newmat,PetscInt sizesset) 1910 { 1911 PetscErrorCode ierr; 1912 PetscMPIInt rank,size; 1913 PetscInt *rowners,i,m,nz,j; 1914 PetscScalar *array,*vals,*vals_ptr; 1915 1916 PetscFunctionBegin; 1917 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1918 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1919 1920 /* determine ownership of all rows */ 1921 if (newmat->rmap->n < 0) m = M/size + ((M % size) > rank); 1922 else m = newmat->rmap->n; 1923 ierr = PetscMalloc((size+2)*sizeof(PetscInt),&rowners);CHKERRQ(ierr); 1924 ierr = MPI_Allgather(&m,1,MPIU_INT,rowners+1,1,MPIU_INT,comm);CHKERRQ(ierr); 1925 rowners[0] = 0; 1926 for (i=2; i<=size; i++) { 1927 rowners[i] += rowners[i-1]; 1928 } 1929 1930 if (!sizesset) { 1931 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 1932 } 1933 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1934 ierr = MatGetArray(newmat,&array);CHKERRQ(ierr); 1935 1936 if (!rank) { 1937 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1938 1939 /* read in my part of the matrix numerical values */ 1940 ierr = PetscBinaryRead(fd,vals,m*N,PETSC_SCALAR);CHKERRQ(ierr); 1941 1942 /* insert into matrix-by row (this is why cannot directly read into array */ 1943 vals_ptr = vals; 1944 for (i=0; i<m; i++) { 1945 for (j=0; j<N; j++) { 1946 array[i + j*m] = *vals_ptr++; 1947 } 1948 } 1949 1950 /* read in other processors and ship out */ 1951 for (i=1; i<size; i++) { 1952 nz = (rowners[i+1] - rowners[i])*N; 1953 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 1954 ierr = MPILong_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1955 } 1956 } else { 1957 /* receive numeric values */ 1958 ierr = PetscMalloc(m*N*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 1959 1960 /* receive message of values*/ 1961 ierr = MPILong_Recv(vals,m*N,MPIU_SCALAR,0,((PetscObject)(newmat))->tag,comm);CHKERRQ(ierr); 1962 1963 /* insert into matrix-by row (this is why cannot directly read into array */ 1964 vals_ptr = vals; 1965 for (i=0; i<m; i++) { 1966 for (j=0; j<N; j++) { 1967 array[i + j*m] = *vals_ptr++; 1968 } 1969 } 1970 } 1971 ierr = PetscFree(rowners);CHKERRQ(ierr); 1972 ierr = PetscFree(vals);CHKERRQ(ierr); 1973 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1974 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1975 PetscFunctionReturn(0); 1976 } 1977 1978 #undef __FUNCT__ 1979 #define __FUNCT__ "MatLoad_MPIDense" 1980 PetscErrorCode MatLoad_MPIDense(Mat newmat,PetscViewer viewer) 1981 { 1982 PetscScalar *vals,*svals; 1983 MPI_Comm comm = ((PetscObject)viewer)->comm; 1984 MPI_Status status; 1985 PetscMPIInt rank,size,tag = ((PetscObject)viewer)->tag,*rowners,*sndcounts,m,maxnz; 1986 PetscInt header[4],*rowlengths = 0,M,N,*cols; 1987 PetscInt *ourlens,*procsnz = 0,*offlens,jj,*mycols,*smycols; 1988 PetscInt i,nz,j,rstart,rend,sizesset=1,grows,gcols; 1989 int fd; 1990 PetscErrorCode ierr; 1991 1992 PetscFunctionBegin; 1993 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1994 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 1995 if (!rank) { 1996 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1997 ierr = PetscBinaryRead(fd,(char *)header,4,PETSC_INT);CHKERRQ(ierr); 1998 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object"); 1999 } 2000 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) sizesset = 0; 2001 2002 ierr = MPI_Bcast(header+1,3,MPIU_INT,0,comm);CHKERRQ(ierr); 2003 M = header[1]; N = header[2]; nz = header[3]; 2004 2005 /* If global rows/cols are set to PETSC_DECIDE, set it to the sizes given in the file */ 2006 if (sizesset && newmat->rmap->N < 0) newmat->rmap->N = M; 2007 if (sizesset && newmat->cmap->N < 0) newmat->cmap->N = N; 2008 2009 /* If global sizes are set, check if they are consistent with that given in the file */ 2010 if (sizesset) { 2011 ierr = MatGetSize(newmat,&grows,&gcols);CHKERRQ(ierr); 2012 } 2013 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); 2014 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); 2015 2016 /* 2017 Handle case where matrix is stored on disk as a dense matrix 2018 */ 2019 if (nz == MATRIX_BINARY_FORMAT_DENSE) { 2020 ierr = MatLoad_MPIDense_DenseInFile(comm,fd,M,N,newmat,sizesset);CHKERRQ(ierr); 2021 PetscFunctionReturn(0); 2022 } 2023 2024 /* determine ownership of all rows */ 2025 if (newmat->rmap->n < 0) m = PetscMPIIntCast(M/size + ((M % size) > rank)); 2026 else m = PetscMPIIntCast(newmat->rmap->n); 2027 ierr = PetscMalloc((size+2)*sizeof(PetscMPIInt),&rowners);CHKERRQ(ierr); 2028 ierr = MPI_Allgather(&m,1,MPI_INT,rowners+1,1,MPI_INT,comm);CHKERRQ(ierr); 2029 rowners[0] = 0; 2030 for (i=2; i<=size; i++) { 2031 rowners[i] += rowners[i-1]; 2032 } 2033 rstart = rowners[rank]; 2034 rend = rowners[rank+1]; 2035 2036 /* distribute row lengths to all processors */ 2037 ierr = PetscMalloc2(rend-rstart,PetscInt,&ourlens,rend-rstart,PetscInt,&offlens);CHKERRQ(ierr); 2038 if (!rank) { 2039 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2040 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2041 ierr = PetscMalloc(size*sizeof(PetscMPIInt),&sndcounts);CHKERRQ(ierr); 2042 for (i=0; i<size; i++) sndcounts[i] = rowners[i+1] - rowners[i]; 2043 ierr = MPI_Scatterv(rowlengths,sndcounts,rowners,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2044 ierr = PetscFree(sndcounts);CHKERRQ(ierr); 2045 } else { 2046 ierr = MPI_Scatterv(0,0,0,MPIU_INT,ourlens,rend-rstart,MPIU_INT,0,comm);CHKERRQ(ierr); 2047 } 2048 2049 if (!rank) { 2050 /* calculate the number of nonzeros on each processor */ 2051 ierr = PetscMalloc(size*sizeof(PetscInt),&procsnz);CHKERRQ(ierr); 2052 ierr = PetscMemzero(procsnz,size*sizeof(PetscInt));CHKERRQ(ierr); 2053 for (i=0; i<size; i++) { 2054 for (j=rowners[i]; j< rowners[i+1]; j++) { 2055 procsnz[i] += rowlengths[j]; 2056 } 2057 } 2058 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2059 2060 /* determine max buffer needed and allocate it */ 2061 maxnz = 0; 2062 for (i=0; i<size; i++) { 2063 maxnz = PetscMax(maxnz,procsnz[i]); 2064 } 2065 ierr = PetscMalloc(maxnz*sizeof(PetscInt),&cols);CHKERRQ(ierr); 2066 2067 /* read in my part of the matrix column indices */ 2068 nz = procsnz[0]; 2069 ierr = PetscMalloc(nz*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2070 ierr = PetscBinaryRead(fd,mycols,nz,PETSC_INT);CHKERRQ(ierr); 2071 2072 /* read in every one elses and ship off */ 2073 for (i=1; i<size; i++) { 2074 nz = procsnz[i]; 2075 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 2076 ierr = MPI_Send(cols,nz,MPIU_INT,i,tag,comm);CHKERRQ(ierr); 2077 } 2078 ierr = PetscFree(cols);CHKERRQ(ierr); 2079 } else { 2080 /* determine buffer space needed for message */ 2081 nz = 0; 2082 for (i=0; i<m; i++) { 2083 nz += ourlens[i]; 2084 } 2085 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&mycols);CHKERRQ(ierr); 2086 2087 /* receive message of column indices*/ 2088 ierr = MPI_Recv(mycols,nz,MPIU_INT,0,tag,comm,&status);CHKERRQ(ierr); 2089 ierr = MPI_Get_count(&status,MPIU_INT,&maxnz);CHKERRQ(ierr); 2090 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2091 } 2092 2093 /* loop over local rows, determining number of off diagonal entries */ 2094 ierr = PetscMemzero(offlens,m*sizeof(PetscInt));CHKERRQ(ierr); 2095 jj = 0; 2096 for (i=0; i<m; i++) { 2097 for (j=0; j<ourlens[i]; j++) { 2098 if (mycols[jj] < rstart || mycols[jj] >= rend) offlens[i]++; 2099 jj++; 2100 } 2101 } 2102 2103 /* create our matrix */ 2104 for (i=0; i<m; i++) { 2105 ourlens[i] -= offlens[i]; 2106 } 2107 2108 if (!sizesset) { 2109 ierr = MatSetSizes(newmat,m,PETSC_DECIDE,M,N);CHKERRQ(ierr); 2110 } 2111 ierr = MatMPIDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 2112 for (i=0; i<m; i++) { 2113 ourlens[i] += offlens[i]; 2114 } 2115 2116 if (!rank) { 2117 ierr = PetscMalloc(maxnz*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2118 2119 /* read in my part of the matrix numerical values */ 2120 nz = procsnz[0]; 2121 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2122 2123 /* insert into matrix */ 2124 jj = rstart; 2125 smycols = mycols; 2126 svals = vals; 2127 for (i=0; i<m; i++) { 2128 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2129 smycols += ourlens[i]; 2130 svals += ourlens[i]; 2131 jj++; 2132 } 2133 2134 /* read in other processors and ship out */ 2135 for (i=1; i<size; i++) { 2136 nz = procsnz[i]; 2137 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 2138 ierr = MPI_Send(vals,nz,MPIU_SCALAR,i,((PetscObject)newmat)->tag,comm);CHKERRQ(ierr); 2139 } 2140 ierr = PetscFree(procsnz);CHKERRQ(ierr); 2141 } else { 2142 /* receive numeric values */ 2143 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 2144 2145 /* receive message of values*/ 2146 ierr = MPI_Recv(vals,nz,MPIU_SCALAR,0,((PetscObject)newmat)->tag,comm,&status);CHKERRQ(ierr); 2147 ierr = MPI_Get_count(&status,MPIU_SCALAR,&maxnz);CHKERRQ(ierr); 2148 if (maxnz != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"something is wrong with file"); 2149 2150 /* insert into matrix */ 2151 jj = rstart; 2152 smycols = mycols; 2153 svals = vals; 2154 for (i=0; i<m; i++) { 2155 ierr = MatSetValues(newmat,1,&jj,ourlens[i],smycols,svals,INSERT_VALUES);CHKERRQ(ierr); 2156 smycols += ourlens[i]; 2157 svals += ourlens[i]; 2158 jj++; 2159 } 2160 } 2161 ierr = PetscFree2(ourlens,offlens);CHKERRQ(ierr); 2162 ierr = PetscFree(vals);CHKERRQ(ierr); 2163 ierr = PetscFree(mycols);CHKERRQ(ierr); 2164 ierr = PetscFree(rowners);CHKERRQ(ierr); 2165 2166 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2167 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2168 PetscFunctionReturn(0); 2169 } 2170 2171 #undef __FUNCT__ 2172 #define __FUNCT__ "MatEqual_MPIDense" 2173 PetscErrorCode MatEqual_MPIDense(Mat A,Mat B,PetscBool *flag) 2174 { 2175 Mat_MPIDense *matB = (Mat_MPIDense*)B->data,*matA = (Mat_MPIDense*)A->data; 2176 Mat a,b; 2177 PetscBool flg; 2178 PetscErrorCode ierr; 2179 2180 PetscFunctionBegin; 2181 a = matA->A; 2182 b = matB->A; 2183 ierr = MatEqual(a,b,&flg);CHKERRQ(ierr); 2184 ierr = MPI_Allreduce(&flg,flag,1,MPI_INT,MPI_LAND,((PetscObject)A)->comm);CHKERRQ(ierr); 2185 PetscFunctionReturn(0); 2186 } 2187 2188 #if defined(PETSC_HAVE_PLAPACK) 2189 2190 #undef __FUNCT__ 2191 #define __FUNCT__ "PetscPLAPACKFinalizePackage" 2192 /*@C 2193 PetscPLAPACKFinalizePackage - This function destroys everything in the Petsc interface to PLAPACK. 2194 Level: developer 2195 2196 .keywords: Petsc, destroy, package, PLAPACK 2197 .seealso: PetscFinalize() 2198 @*/ 2199 PetscErrorCode PetscPLAPACKFinalizePackage(void) 2200 { 2201 PetscErrorCode ierr; 2202 2203 PetscFunctionBegin; 2204 ierr = PLA_Finalize();CHKERRQ(ierr); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 #undef __FUNCT__ 2209 #define __FUNCT__ "PetscPLAPACKInitializePackage" 2210 /*@C 2211 PetscPLAPACKInitializePackage - This function initializes everything in the Petsc interface to PLAPACK. It is 2212 called from MatCreate_MPIDense() the first time an MPI dense matrix is called. 2213 2214 Input Parameter: 2215 . comm - the communicator the matrix lives on 2216 2217 Level: developer 2218 2219 Notes: PLAPACK does not have a good fit with MPI communicators; all (parallel) PLAPACK objects have to live in the 2220 same communicator (because there is some global state that is initialized and used for all matrices). In addition if 2221 PLAPACK is initialized (that is the initial matrices created) are on subcommunicators of MPI_COMM_WORLD, these subcommunicators 2222 cannot overlap. 2223 2224 .keywords: Petsc, initialize, package, PLAPACK 2225 .seealso: PetscSysInitializePackage(), PetscInitialize() 2226 @*/ 2227 PetscErrorCode PetscPLAPACKInitializePackage(MPI_Comm comm) 2228 { 2229 PetscMPIInt size; 2230 PetscErrorCode ierr; 2231 2232 PetscFunctionBegin; 2233 if (!PLA_Initialized(PETSC_NULL)) { 2234 2235 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2236 Plapack_nprows = 1; 2237 Plapack_npcols = size; 2238 2239 ierr = PetscOptionsBegin(comm,PETSC_NULL,"PLAPACK Options","Mat");CHKERRQ(ierr); 2240 ierr = PetscOptionsInt("-mat_plapack_nprows","row dimension of 2D processor mesh","None",Plapack_nprows,&Plapack_nprows,PETSC_NULL);CHKERRQ(ierr); 2241 ierr = PetscOptionsInt("-mat_plapack_npcols","column dimension of 2D processor mesh","None",Plapack_npcols,&Plapack_npcols,PETSC_NULL);CHKERRQ(ierr); 2242 #if defined(PETSC_USE_DEBUG) 2243 Plapack_ierror = 3; 2244 #else 2245 Plapack_ierror = 0; 2246 #endif 2247 ierr = PetscOptionsInt("-mat_plapack_ckerror","error checking flag","None",Plapack_ierror,&Plapack_ierror,PETSC_NULL);CHKERRQ(ierr); 2248 if (Plapack_ierror){ 2249 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_TRUE,PETSC_TRUE,PETSC_FALSE );CHKERRQ(ierr); 2250 } else { 2251 ierr = PLA_Set_error_checking(Plapack_ierror,PETSC_FALSE,PETSC_FALSE,PETSC_FALSE );CHKERRQ(ierr); 2252 } 2253 2254 Plapack_nb_alg = 0; 2255 ierr = PetscOptionsInt("-mat_plapack_nb_alg","algorithmic block size","None",Plapack_nb_alg,&Plapack_nb_alg,PETSC_NULL);CHKERRQ(ierr); 2256 if (Plapack_nb_alg) { 2257 ierr = pla_Environ_set_nb_alg (PLA_OP_ALL_ALG,Plapack_nb_alg);CHKERRQ(ierr); 2258 } 2259 PetscOptionsEnd(); 2260 2261 ierr = PLA_Comm_1D_to_2D(comm,Plapack_nprows,Plapack_npcols,&Plapack_comm_2d);CHKERRQ(ierr); 2262 ierr = PLA_Init(Plapack_comm_2d);CHKERRQ(ierr); 2263 ierr = PetscRegisterFinalize(PetscPLAPACKFinalizePackage);CHKERRQ(ierr); 2264 } 2265 PetscFunctionReturn(0); 2266 } 2267 2268 #endif 2269