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