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