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