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