本文整理汇总了C++中MatCreate函数的典型用法代码示例。如果您正苦于以下问题:C++ MatCreate函数的具体用法?C++ MatCreate怎么用?C++ MatCreate使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MatCreate函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc,char **argv)
{
Mat A,A11,A12,A21,A22;
Vec X,X1,X2,Y,Z,Z1,Z2;
PetscScalar *a,*b,*x,*y,*z,v,one=1;
PetscReal nrm;
PetscErrorCode ierr;
PetscInt size=8,size1=6,size2=2, i,j;
PetscInitialize(&argc,&argv,0,help);
/*
* Create matrix and three vectors: these are all normal
*/
ierr = PetscMalloc(size*size*sizeof(PetscScalar),&a);CHKERRQ(ierr);
ierr = PetscMalloc(size*size*sizeof(PetscScalar),&b);CHKERRQ(ierr);
for (i=0; i<size; i++) {
for (j=0; j<size; j++) {
a[i+j*size] = rand(); b[i+j*size] = a[i+j*size];
}
}
ierr = MatCreate(MPI_COMM_SELF,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,size,size,size,size);CHKERRQ(ierr);
ierr = MatSetType(A,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(A,a);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = PetscMalloc(size*sizeof(PetscScalar),&x);CHKERRQ(ierr);
for (i=0; i<size; i++) {
x[i] = one;
}
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,x,&X);CHKERRQ(ierr);
ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
ierr = VecAssemblyEnd(X);CHKERRQ(ierr);
ierr = PetscMalloc(size*sizeof(PetscScalar),&y);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,y,&Y);CHKERRQ(ierr);
ierr = VecAssemblyBegin(Y);CHKERRQ(ierr);
ierr = VecAssemblyEnd(Y);CHKERRQ(ierr);
ierr = PetscMalloc(size*sizeof(PetscScalar),&z);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size,z,&Z);CHKERRQ(ierr);
ierr = VecAssemblyBegin(Z);CHKERRQ(ierr);
ierr = VecAssemblyEnd(Z);CHKERRQ(ierr);
/*
* Now create submatrices and subvectors
*/
ierr = MatCreate(MPI_COMM_SELF,&A11);CHKERRQ(ierr);
ierr = MatSetSizes(A11,size1,size1,size1,size1);CHKERRQ(ierr);
ierr = MatSetType(A11,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(A11,b);CHKERRQ(ierr);
ierr = MatSeqDenseSetLDA(A11,size);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A11,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A11,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatCreate(MPI_COMM_SELF,&A12);CHKERRQ(ierr);
ierr = MatSetSizes(A12,size1,size2,size1,size2);CHKERRQ(ierr);
ierr = MatSetType(A12,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(A12,b+size1*size);CHKERRQ(ierr);
ierr = MatSeqDenseSetLDA(A12,size);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A12,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A12,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatCreate(MPI_COMM_SELF,&A21);CHKERRQ(ierr);
ierr = MatSetSizes(A21,size2,size1,size2,size1);CHKERRQ(ierr);
ierr = MatSetType(A21,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(A21,b+size1);CHKERRQ(ierr);
ierr = MatSeqDenseSetLDA(A21,size);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A21,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A21,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatCreate(MPI_COMM_SELF,&A22);CHKERRQ(ierr);
ierr = MatSetSizes(A22,size2,size2,size2,size2);CHKERRQ(ierr);
ierr = MatSetType(A22,MATSEQDENSE);CHKERRQ(ierr);
ierr = MatSeqDenseSetPreallocation(A22,b+size1*size+size1);CHKERRQ(ierr);
ierr = MatSeqDenseSetLDA(A22,size);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A22,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A22,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,x,&X1);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,x+size1,&X2);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size1,z,&Z1);CHKERRQ(ierr);
ierr = VecCreateSeqWithArray(MPI_COMM_SELF,1,size2,z+size1,&Z2);CHKERRQ(ierr);
/*
* Now multiple matrix times input in two ways;
* compare the results
*/
ierr = MatMult(A,X,Y);CHKERRQ(ierr);
ierr = MatMult(A11,X1,Z1);CHKERRQ(ierr);
ierr = MatMultAdd(A12,X2,Z1,Z1);CHKERRQ(ierr);
ierr = MatMult(A22,X2,Z2);CHKERRQ(ierr);
ierr = MatMultAdd(A21,X1,Z2,Z2);CHKERRQ(ierr);
ierr = VecAXPY(Z,-1.0,Y);CHKERRQ(ierr);
ierr = VecNorm(Z,NORM_2,&nrm);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Test1; error norm=%G\n",nrm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"MatMult the usual way:\n");CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:feelpp,项目名称:debian-petsc,代码行数:101,代码来源:ex2.c
示例2: main
int main(int argc,char **argv)
{
Mat M,C,K,A[3]; /* problem matrices */
PEP pep; /* polynomial eigenproblem solver context */
PetscInt m=6,n,II,Istart,Iend,i,j;
PetscScalar z=1.0;
PetscReal h;
char str[50];
PetscErrorCode ierr;
SlepcInitialize(&argc,&argv,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
if (m<2) SETERRQ(PETSC_COMM_SELF,1,"m must be at least 2");
ierr = PetscOptionsGetScalar(NULL,"-z",&z,NULL);CHKERRQ(ierr);
h = 1.0/m;
n = m*(m-1);
ierr = SlepcSNPrintfScalar(str,50,z,PETSC_FALSE);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAcoustic wave 2-D, n=%D (m=%D), z=%s\n\n",n,m,str);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* K has a pattern similar to the 2D Laplacian */
ierr = MatCreate(PETSC_COMM_WORLD,&K);CHKERRQ(ierr);
ierr = MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(K);CHKERRQ(ierr);
ierr = MatSetUp(K);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(K,&Istart,&Iend);CHKERRQ(ierr);
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
if (i>0) { ierr = MatSetValue(K,II,II-m,(j==m-1)?-0.5:-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-2) { ierr = MatSetValue(K,II,II+m,(j==m-1)?-0.5:-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if (j>0) { ierr = MatSetValue(K,II,II-1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(K,II,II+1,-1.0,INSERT_VALUES);CHKERRQ(ierr); }
ierr = MatSetValue(K,II,II,(j==m-1)?2.0:4.0,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* C is the zero matrix except for a few nonzero elements on the diagonal */
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(C);CHKERRQ(ierr);
ierr = MatSetUp(C);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(C,&Istart,&Iend);CHKERRQ(ierr);
for (i=Istart;i<Iend;i++) {
if (i%m==m-1) {
ierr = MatSetValue(C,i,i,-2*PETSC_PI*h/z,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* M is a diagonal matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&M);CHKERRQ(ierr);
ierr = MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(M);CHKERRQ(ierr);
ierr = MatSetUp(M);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(M,&Istart,&Iend);CHKERRQ(ierr);
for (i=Istart;i<Iend;i++) {
if (i%m==m-1) {
ierr = MatSetValue(M,i,i,2*PETSC_PI*PETSC_PI*h*h,INSERT_VALUES);CHKERRQ(ierr);
} else {
ierr = MatSetValue(M,i,i,4*PETSC_PI*PETSC_PI*h*h,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and solve the problem
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = PEPCreate(PETSC_COMM_WORLD,&pep);CHKERRQ(ierr);
A[0] = K; A[1] = C; A[2] = M;
ierr = PEPSetOperators(pep,3,A);CHKERRQ(ierr);
ierr = PEPSetFromOptions(pep);CHKERRQ(ierr);
ierr = PEPSolve(pep);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = PEPPrintSolution(pep,NULL);CHKERRQ(ierr);
ierr = PEPDestroy(&pep);CHKERRQ(ierr);
ierr = MatDestroy(&M);CHKERRQ(ierr);
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = MatDestroy(&K);CHKERRQ(ierr);
ierr = SlepcFinalize();CHKERRQ(ierr);
return 0;
}
开发者ID:OpenCMISS-Dependencies,项目名称:slepc,代码行数:97,代码来源:acoustic_wave_2d.c
示例3: main
int main(int argc,char **argv)
{
Mat A[NMAT]; /* problem matrices */
PEP pep; /* polynomial eigenproblem solver context */
PetscInt n,m=8,k,II,Istart,Iend,i,j;
PetscReal c[10] = { 0.6, 1.3, 1.3, 0.1, 0.1, 1.2, 1.0, 1.0, 1.2, 1.0 };
PetscBool flg;
PetscErrorCode ierr;
SlepcInitialize(&argc,&argv,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
n = m*m;
k = 10;
ierr = PetscOptionsGetRealArray(NULL,"-c",c,&k,&flg);CHKERRQ(ierr);
if (flg && k!=10) SETERRQ1(PETSC_COMM_WORLD,1,"The number of parameters -c should be 10, you provided %D",k);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nButterfly problem, n=%D (m=%D)\n\n",n,m);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the polynomial matrices
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* initialize matrices */
for (i=0;i<NMAT;i++) {
ierr = MatCreate(PETSC_COMM_WORLD,&A[i]);CHKERRQ(ierr);
ierr = MatSetSizes(A[i],PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A[i]);CHKERRQ(ierr);
ierr = MatSetUp(A[i]);CHKERRQ(ierr);
}
ierr = MatGetOwnershipRange(A[0],&Istart,&Iend);CHKERRQ(ierr);
/* A0 */
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
ierr = MatSetValue(A[0],II,II,4.0*c[0]/6.0+4.0*c[1]/6.0,INSERT_VALUES);CHKERRQ(ierr);
if (j>0) { ierr = MatSetValue(A[0],II,II-1,c[0]/6.0,INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(A[0],II,II+1,c[0]/6.0,INSERT_VALUES);CHKERRQ(ierr); }
if (i>0) { ierr = MatSetValue(A[0],II,II-m,c[1]/6.0,INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-1) { ierr = MatSetValue(A[0],II,II+m,c[1]/6.0,INSERT_VALUES);CHKERRQ(ierr); }
}
/* A1 */
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
if (j>0) { ierr = MatSetValue(A[1],II,II-1,c[2],INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(A[1],II,II+1,-c[2],INSERT_VALUES);CHKERRQ(ierr); }
if (i>0) { ierr = MatSetValue(A[1],II,II-m,c[3],INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-1) { ierr = MatSetValue(A[1],II,II+m,-c[3],INSERT_VALUES);CHKERRQ(ierr); }
}
/* A2 */
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
ierr = MatSetValue(A[2],II,II,-2.0*c[4]-2.0*c[5],INSERT_VALUES);CHKERRQ(ierr);
if (j>0) { ierr = MatSetValue(A[2],II,II-1,c[4],INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(A[2],II,II+1,c[4],INSERT_VALUES);CHKERRQ(ierr); }
if (i>0) { ierr = MatSetValue(A[2],II,II-m,c[5],INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-1) { ierr = MatSetValue(A[2],II,II+m,c[5],INSERT_VALUES);CHKERRQ(ierr); }
}
/* A3 */
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
if (j>0) { ierr = MatSetValue(A[3],II,II-1,c[6],INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(A[3],II,II+1,-c[6],INSERT_VALUES);CHKERRQ(ierr); }
if (i>0) { ierr = MatSetValue(A[3],II,II-m,c[7],INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-1) { ierr = MatSetValue(A[3],II,II+m,-c[7],INSERT_VALUES);CHKERRQ(ierr); }
}
/* A4 */
for (II=Istart;II<Iend;II++) {
i = II/m; j = II-i*m;
ierr = MatSetValue(A[4],II,II,2.0*c[8]+2.0*c[9],INSERT_VALUES);CHKERRQ(ierr);
if (j>0) { ierr = MatSetValue(A[4],II,II-1,-c[8],INSERT_VALUES);CHKERRQ(ierr); }
if (j<m-1) { ierr = MatSetValue(A[4],II,II+1,-c[8],INSERT_VALUES);CHKERRQ(ierr); }
if (i>0) { ierr = MatSetValue(A[4],II,II-m,-c[9],INSERT_VALUES);CHKERRQ(ierr); }
if (i<m-1) { ierr = MatSetValue(A[4],II,II+m,-c[9],INSERT_VALUES);CHKERRQ(ierr); }
}
/* assemble matrices */
for (i=0;i<NMAT;i++) {
ierr = MatAssemblyBegin(A[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
for (i=0;i<NMAT;i++) {
ierr = MatAssemblyEnd(A[i],MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the eigensolver and solve the problem
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = PEPCreate(PETSC_COMM_WORLD,&pep);CHKERRQ(ierr);
ierr = PEPSetOperators(pep,NMAT,A);CHKERRQ(ierr);
ierr = PEPSetFromOptions(pep);CHKERRQ(ierr);
ierr = PEPSolve(pep);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//.........这里部分代码省略.........
开发者ID:OpenCMISS-Dependencies,项目名称:slepc,代码行数:101,代码来源:butterfly.c
示例4: main
int main(int argc,char **argv)
{
Mat A,B,C,D;
PetscInt i,M=10,N=5,j,nrows,ncols,am,an,rstart,rend;
PetscErrorCode ierr;
PetscRandom r;
PetscBool equal,iselemental;
PetscReal fill = 1.0;
IS isrows,iscols;
const PetscInt *rows,*cols;
PetscScalar *v,rval;
#if defined(PETSC_HAVE_ELEMENTAL)
PetscBool Test_MatMatMult=PETSC_TRUE;
#else
PetscBool Test_MatMatMult=PETSC_FALSE;
#endif
PetscMPIInt size;
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr);
ierr = MatSetType(A,MATDENSE);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&r);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
/* Set local matrix entries */
ierr = MatGetOwnershipIS(A,&isrows,&iscols);CHKERRQ(ierr);
ierr = ISGetLocalSize(isrows,&nrows);CHKERRQ(ierr);
ierr = ISGetIndices(isrows,&rows);CHKERRQ(ierr);
ierr = ISGetLocalSize(iscols,&ncols);CHKERRQ(ierr);
ierr = ISGetIndices(iscols,&cols);CHKERRQ(ierr);
ierr = PetscMalloc1(nrows*ncols,&v);CHKERRQ(ierr);
for (i=0; i<nrows; i++) {
for (j=0; j<ncols; j++) {
ierr = PetscRandomGetValue(r,&rval);CHKERRQ(ierr);
v[i*ncols+j] = rval;
}
}
ierr = MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = ISRestoreIndices(isrows,&rows);CHKERRQ(ierr);
ierr = ISRestoreIndices(iscols,&cols);CHKERRQ(ierr);
ierr = ISDestroy(&isrows);CHKERRQ(ierr);
ierr = ISDestroy(&iscols);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
/* Test MatTranspose() */
ierr = MatCreateTranspose(A,&C);CHKERRQ(ierr);
ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); /* B = A^T */
ierr = MatMultEqual(C,B,10,&equal);CHKERRQ(ierr);
if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
ierr = MatTranspose(A,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr); /* B = A^T */
ierr = MatMultEqual(C,B,10,&equal);CHKERRQ(ierr);
if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
ierr = MatDestroy(&B);CHKERRQ(ierr);
ierr = MatDuplicate(A,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
ierr = MatTranspose(B,MAT_INPLACE_MATRIX,&B);CHKERRQ(ierr);
ierr = MatMultEqual(C,B,10,&equal);CHKERRQ(ierr);
if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
ierr = MatDestroy(&B);CHKERRQ(ierr);
ierr = MatDestroy(&C);CHKERRQ(ierr);
/* Test MatMatMult() */
if (Test_MatMatMult) {
#if !defined(PETSC_HAVE_ELEMENTAL)
if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL");
#endif
ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); /* B = A^T */
ierr = MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); /* C = B*A = A^T*A */
ierr = MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr);
/* Test MatDuplicate for matrix product */
ierr = MatDuplicate(C,MAT_COPY_VALUES,&D);CHKERRQ(ierr);
ierr = MatDestroy(&D);CHKERRQ(ierr);
/* Test B*A*x = C*x for n random vector x */
ierr = MatMatMultEqual(B,A,C,10,&equal);CHKERRQ(ierr);
if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = MatMatMultSymbolic(B,A,fill,&C);CHKERRQ(ierr);
for (i=0; i<2; i++) {
/* Repeat the numeric product to test reuse of the previous symbolic product */
ierr = MatMatMultNumeric(B,A,C);CHKERRQ(ierr);
ierr = MatMatMultEqual(B,A,C,10,&equal);CHKERRQ(ierr);
if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
}
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = MatDestroy(&B);CHKERRQ(ierr);
}
/* Test MatTransposeMatMult() */
//.........这里部分代码省略.........
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:ex104.c
示例5: main
int main(int argc,char **args)
{
Mat C;
Vec u,b;
PetscErrorCode ierr;
PetscMPIInt size,rank;
PetscInt i,m = 5,N,start,end,M,idx[4];
PetscInt j,nrsub,ncsub,*rsub,*csub,mystart,myend;
PetscBool flg;
PetscScalar one = 1.0,Ke[16],*vals;
PetscReal h,norm;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);
N = (m+1)*(m+1); /* dimension of matrix */
M = m*m; /* number of elements */
h = 1.0/m; /* mesh width */
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
/* Create stiffness matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetFromOptions(C);CHKERRQ(ierr);
ierr = MatSetUp(C);CHKERRQ(ierr);
start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
end = start + M/size + ((M%size) > rank);
/* Form the element stiffness for the Laplacian */
ierr = FormElementStiffness(h*h,Ke);CHKERRQ(ierr);
for (i=start; i<end; i++) {
/* location of lower left corner of element */
/* node numbers for the four corners of element */
idx[0] = (m+1)*(i/m) + (i % m);
idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
ierr = MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Assemble the matrix again */
ierr = MatZeroEntries(C);CHKERRQ(ierr);
for (i=start; i<end; i++) {
/* location of lower left corner of element */
/* node numbers for the four corners of element */
idx[0] = (m+1)*(i/m) + (i % m);
idx[1] = idx[0]+1; idx[2] = idx[1] + m + 1; idx[3] = idx[2] - 1;
ierr = MatSetValues(C,4,idx,4,idx,Ke,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Create test vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecSet(u,one);CHKERRQ(ierr);
/* Check error */
ierr = MatMult(C,u,b);CHKERRQ(ierr);
ierr = VecNorm(b,NORM_2,&norm);CHKERRQ(ierr);
if (norm > PETSC_SQRT_MACHINE_EPSILON) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error b %g should be near 0\n",(double)norm);CHKERRQ(ierr);
}
/* Now test MatGetValues() */
ierr = PetscOptionsHasName(NULL,NULL,"-get_values",&flg);CHKERRQ(ierr);
if (flg) {
ierr = MatGetOwnershipRange(C,&mystart,&myend);CHKERRQ(ierr);
nrsub = myend - mystart; ncsub = 4;
ierr = PetscMalloc1(nrsub*ncsub,&vals);CHKERRQ(ierr);
ierr = PetscMalloc1(nrsub,&rsub);CHKERRQ(ierr);
ierr = PetscMalloc1(ncsub,&csub);CHKERRQ(ierr);
for (i=myend-1; i>=mystart; i--) rsub[myend-i-1] = i;
for (i=0; i<ncsub; i++) csub[i] = 2*(ncsub-i) + mystart;
ierr = MatGetValues(C,nrsub,rsub,ncsub,csub,vals);CHKERRQ(ierr);
ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"processor number %d: start=%D, end=%D, mystart=%D, myend=%D\n",rank,start,end,mystart,myend);CHKERRQ(ierr);
for (i=0; i<nrsub; i++) {
for (j=0; j<ncsub; j++) {
if (PetscImaginaryPart(vals[i*ncsub+j]) != 0.0) {
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD," C[%D, %D] = %g + %g i\n",rsub[i],csub[j],(double)PetscRealPart(vals[i*ncsub+j]),(double)PetscImaginaryPart(vals[i*ncsub+j]));CHKERRQ(ierr);
} else {
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD," C[%D, %D] = %g\n",rsub[i],csub[j],(double)PetscRealPart(vals[i*ncsub+j]));CHKERRQ(ierr);
}
}
}
ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
ierr = PetscFree(rsub);CHKERRQ(ierr);
ierr = PetscFree(csub);CHKERRQ(ierr);
ierr = PetscFree(vals);CHKERRQ(ierr);
}
/* Free data structures */
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:ex19.c
示例6: main
int main(int argc,char **argv)
{
TS ts; /* nonlinear solver */
Vec x; /* solution, residual vectors */
Mat A; /* Jacobian matrix */
Mat Jacp; /* JacobianP matrix */
PetscInt steps;
PetscReal ftime =0.5;
PetscBool monitor = PETSC_FALSE;
PetscScalar *x_ptr;
PetscMPIInt size;
struct _n_User user;
PetscErrorCode ierr;
Vec lambda[2],mu[2];
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Initialize program
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
PetscInitialize(&argc,&argv,NULL,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
user.mu = 1;
user.next_output = 0.0;
ierr = PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create necessary matrix and vectors, solve same ODE on every process
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&Jacp);CHKERRQ(ierr);
ierr = MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr);
ierr = MatSetFromOptions(Jacp);CHKERRQ(ierr);
ierr = MatSetUp(Jacp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create timestepping solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
ierr = TSSetType(ts,TSRK);CHKERRQ(ierr);
ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&user);CHKERRQ(ierr);
ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
if (monitor) {
ierr = TSMonitorSet(ts,Monitor,&user,NULL);CHKERRQ(ierr);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set initial conditions
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = VecGetArray(x,&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 2; x_ptr[1] = 0.66666654321;
ierr = VecRestoreArray(x,&x_ptr);CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,0.0,.001);CHKERRQ(ierr);
/*
Have the TS save its trajectory so that TSAdjointSolve() may be used
*/
ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSolve(ts,x);CHKERRQ(ierr);
ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,steps,(double)ftime);CHKERRQ(ierr);
ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Start the Adjoint model
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreateVecs(A,&lambda[0],NULL);CHKERRQ(ierr);
ierr = MatCreateVecs(A,&lambda[1],NULL);CHKERRQ(ierr);
/* Reset initial conditions for the adjoint integration */
ierr = VecGetArray(lambda[0],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 1.0; x_ptr[1] = 0.0;
ierr = VecRestoreArray(lambda[0],&x_ptr);CHKERRQ(ierr);
ierr = VecGetArray(lambda[1],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = 0.0; x_ptr[1] = 1.0;
ierr = VecRestoreArray(lambda[1],&x_ptr);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:masa-ito,项目名称:PETScToPoisson,代码行数:101,代码来源:ex16adj.c
示例7: MatPtAPSymbolic_MPIAIJ_MPIAIJ
//.........这里部分代码省略.........
/* add non-zero cols of AP into the sorted linked list lnk */
ierr = PetscLLCondensedAddSorted(apnz,Jptr,lnk,lnkbt);CHKERRQ(ierr);
}
/* add received col data into lnk */
for (k=0; k<merge->nrecv; k++) { /* k-th received message */
if (i == *nextrow[k]) { /* i-th row */
nzi = *(nextci[k]+1) - *nextci[k];
Jptr = buf_rj[k] + *nextci[k];
ierr = PetscLLCondensedAddSorted(nzi,Jptr,lnk,lnkbt);CHKERRQ(ierr);
nextrow[k]++; nextci[k]++;
}
}
nnz = lnk[0];
/* if free space is not available, make more free space */
if (current_space->local_remaining<nnz) {
ierr = PetscFreeSpaceGet(nnz+current_space->total_array_size,¤t_space);CHKERRQ(ierr);
nspacedouble++;
}
/* copy data into free space, then initialize lnk */
ierr = PetscLLCondensedClean(pN,nnz,current_space->array,lnk,lnkbt);CHKERRQ(ierr);
ierr = MatPreallocateSet(i+owners[rank],nnz,current_space->array,dnz,onz);CHKERRQ(ierr);
current_space->array += nnz;
current_space->local_used += nnz;
current_space->local_remaining -= nnz;
pti[i+1] = pti[i] + nnz;
if (nnz > rmax) rmax = nnz;
}
ierr = MatRestoreSymbolicTranspose_SeqAIJ(p->A,&pdti,&pdtj);CHKERRQ(ierr);
ierr = PetscFree3(buf_ri_k,nextrow,nextci);CHKERRQ(ierr);
ierr = PetscMalloc1((pti[pn]+1),&ptj);CHKERRQ(ierr);
ierr = PetscFreeSpaceContiguous(&free_space,ptj);CHKERRQ(ierr);
afill_tmp = (PetscReal)pti[pn]/(pi_loc[pm] + api[am]+1);
if (afill_tmp > afill) afill = afill_tmp;
ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
/* create symbolic parallel matrix Cmpi */
/*--------------------------------------*/
ierr = MatCreate(comm,&Cmpi);CHKERRQ(ierr);
ierr = MatSetSizes(Cmpi,pn,pn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSetBlockSizes(Cmpi,P->cmap->bs,P->cmap->bs);CHKERRQ(ierr);
ierr = MatSetType(Cmpi,MATMPIAIJ);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(Cmpi,0,dnz,0,onz);CHKERRQ(ierr);
ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr);
merge->bi = pti; /* Cseq->i */
merge->bj = ptj; /* Cseq->j */
merge->coi = coi; /* Co->i */
merge->coj = coj; /* Co->j */
merge->buf_ri = buf_ri;
merge->buf_rj = buf_rj;
merge->owners_co = owners_co;
merge->destroy = Cmpi->ops->destroy;
merge->duplicate = Cmpi->ops->duplicate;
/* Cmpi is not ready for use - assembly will be done by MatPtAPNumeric() */
Cmpi->assembled = PETSC_FALSE;
Cmpi->ops->destroy = MatDestroy_MPIAIJ_PtAP;
Cmpi->ops->duplicate = MatDuplicate_MPIAIJ_MatPtAP;
/* attach the supporting struct to Cmpi for reuse */
c = (Mat_MPIAIJ*)Cmpi->data;
c->ptap = ptap;
ptap->api = api;
ptap->apj = apj;
ptap->rmax = ap_rmax;
*C = Cmpi;
/* flag 'scalable' determines which implementations to be used:
0: do dense axpy in MatPtAPNumeric() - fast, but requires storage of a nonscalable dense array apa;
1: do sparse axpy in MatPtAPNumeric() - might slow, uses a sparse array apa */
/* set default scalable */
ptap->scalable = PETSC_TRUE;
ierr = PetscOptionsGetBool(((PetscObject)Cmpi)->prefix,"-matptap_scalable",&ptap->scalable,NULL);CHKERRQ(ierr);
if (!ptap->scalable) { /* Do dense axpy */
ierr = PetscCalloc1(pN,&ptap->apa);CHKERRQ(ierr);
} else {
ierr = PetscCalloc1(ap_rmax+1,&ptap->apa);CHKERRQ(ierr);
}
#if defined(PTAP_PROFILE)
ierr = PetscTime(&t4);CHKERRQ(ierr);
if (rank==1) PetscPrintf(MPI_COMM_SELF," [%d] PtAPSymbolic %g/P + %g/AP + %g/comm + %g/PtAP = %g\n",rank,t1-t0,t2-t1,t3-t2,t4-t3,t4-t0);CHKERRQ(ierr);
#endif
#if defined(PETSC_USE_INFO)
if (pti[pn] != 0) {
ierr = PetscInfo3(Cmpi,"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
ierr = PetscInfo1(Cmpi,"Use MatPtAP(A,P,MatReuse,%G,&C) for best performance.\n",afill);CHKERRQ(ierr);
} else {
ierr = PetscInfo(Cmpi,"Empty matrix product\n");CHKERRQ(ierr);
}
#endif
PetscFunctionReturn(0);
}
开发者ID:ZJLi2013,项目名称:petsc,代码行数:101,代码来源:mpiptap.c
示例8: MatGetMultiProcBlock_MPIAIJ
/*
Developers Note: This is used directly by some preconditioners, hence is PETSC_EXTERN
*/
PETSC_EXTERN PetscErrorCode MatGetMultiProcBlock_MPIAIJ(Mat mat, MPI_Comm subComm, MatReuse scall,Mat *subMat)
{
PetscErrorCode ierr;
Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data;
Mat_SeqAIJ *aijB = (Mat_SeqAIJ*)aij->B->data;
PetscMPIInt commRank,subCommSize,subCommRank;
PetscMPIInt *commRankMap,subRank,rank,commsize;
PetscInt *garrayCMap,col,i,j,*nnz,newRow,newCol;
PetscFunctionBegin;
ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&commsize);CHKERRQ(ierr);
ierr = MPI_Comm_size(subComm,&subCommSize);CHKERRQ(ierr);
/* create subMat object with the relavent layout */
if (scall == MAT_INITIAL_MATRIX) {
ierr = MatCreate(subComm,subMat);CHKERRQ(ierr);
ierr = MatSetType(*subMat,MATMPIAIJ);CHKERRQ(ierr);
ierr = MatSetSizes(*subMat,mat->rmap->n,mat->cmap->n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
ierr = MatSetBlockSizes(*subMat,mat->rmap->bs,mat->cmap->bs);CHKERRQ(ierr);
/* need to setup rmap and cmap before Preallocation */
ierr = PetscLayoutSetBlockSize((*subMat)->rmap,mat->rmap->bs);CHKERRQ(ierr);
ierr = PetscLayoutSetBlockSize((*subMat)->cmap,mat->cmap->bs);CHKERRQ(ierr);
ierr = PetscLayoutSetUp((*subMat)->rmap);CHKERRQ(ierr);
ierr = PetscLayoutSetUp((*subMat)->cmap);CHKERRQ(ierr);
}
/* create a map of comm_rank from subComm to comm - should commRankMap and garrayCMap be kept for reused? */
ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)mat),&commRank);CHKERRQ(ierr);
ierr = MPI_Comm_rank(subComm,&subCommRank);CHKERRQ(ierr);
ierr = PetscMalloc(subCommSize*sizeof(PetscMPIInt),&commRankMap);CHKERRQ(ierr);
ierr = MPI_Allgather(&commRank,1,MPI_INT,commRankMap,1,MPI_INT,subComm);CHKERRQ(ierr);
/* Traverse garray and identify column indices [of offdiag mat] that
should be discarded. For the ones not discarded, store the newCol+1
value in garrayCMap */
ierr = PetscMalloc(aij->B->cmap->n*sizeof(PetscInt),&garrayCMap);CHKERRQ(ierr);
ierr = PetscMemzero(garrayCMap,aij->B->cmap->n*sizeof(PetscInt));CHKERRQ(ierr);
for (i=0; i<aij->B->cmap->n; i++) {
col = aij->garray[i];
for (subRank=0; subRank<subCommSize; subRank++) {
rank = commRankMap[subRank];
if ((col >= mat->cmap->range[rank]) && (col < mat->cmap->range[rank+1])) {
garrayCMap[i] = (*subMat)->cmap->range[subRank] + col - mat->cmap->range[rank]+1;
break;
}
}
}
if (scall == MAT_INITIAL_MATRIX) {
/* Now compute preallocation for the offdiag mat */
ierr = PetscMalloc(aij->B->rmap->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr);
ierr = PetscMemzero(nnz,aij->B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr);
for (i=0; i<aij->B->rmap->n; i++) {
for (j=aijB->i[i]; j<aijB->i[i+1]; j++) {
if (garrayCMap[aijB->j[j]]) nnz[i]++;
}
}
ierr = MatMPIAIJSetPreallocation(*(subMat),0,NULL,0,nnz);CHKERRQ(ierr);
/* reuse diag block with the new submat */
ierr = MatDestroy(&((Mat_MPIAIJ*)((*subMat)->data))->A);CHKERRQ(ierr);
((Mat_MPIAIJ*)((*subMat)->data))->A = aij->A;
ierr = PetscObjectReference((PetscObject)aij->A);CHKERRQ(ierr);
} else if (((Mat_MPIAIJ*)(*subMat)->data)->A != aij->A) {
PetscObject obj = (PetscObject)((Mat_MPIAIJ*)((*subMat)->data))->A;
ierr = PetscObjectReference((PetscObject)obj);CHKERRQ(ierr);
((Mat_MPIAIJ*)((*subMat)->data))->A = aij->A;
ierr = PetscObjectReference((PetscObject)aij->A);CHKERRQ(ierr);
}
/* Now traverse aij->B and insert values into subMat */
for (i=0; i<aij->B->rmap->n; i++) {
newRow = (*subMat)->rmap->range[subCommRank] + i;
for (j=aijB->i[i]; j<aijB->i[i+1]; j++) {
newCol = garrayCMap[aijB->j[j]];
if (newCol) {
newCol--; /* remove the increment */
ierr = MatSetValues(*subMat,1,&newRow,1,&newCol,(aijB->a+j),INSERT_VALUES);CHKERRQ(ierr);
}
}
}
/* assemble the submat */
ierr = MatAssemblyBegin(*subMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(*subMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* deallocate temporary data */
ierr = PetscFree(commRankMap);CHKERRQ(ierr);
ierr = PetscFree(garrayCMap);CHKERRQ(ierr);
if (scall == MAT_INITIAL_MATRIX) {
ierr = PetscFree(nnz);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:hansec,项目名称:petsc,代码行数:101,代码来源:mpb_aij.c
示例9: main
int main(int argc,char **args)
{
/*PETSc Mat Object */
Mat pMat;
/* Input matrix market file and output PETSc binary file */
char inputFile[128],outputFile[128],buf[128];
/* number rows, columns, non zeros etc */
int i,j,m,n,nnz,ierr,col,row;
/*We compute no of nozeros per row for PETSc Mat object pre-allocation*/
int *nnzPtr;
/*Maximum nonzero in nay row */
int maxNNZperRow=0;
/*Row number containing max non zero elements */
int maxRowNum = 0;
/*Just no of comments that will be ignore during successive read of file */
int numComments=0;
PetscScalar zero=0;
/* This is variable of type double */
PetscScalar val;
/*File handle for read and write*/
FILE* file;
/*File handle for writing nonzero elements distribution per row */
FILE *fileRowDist;
/*PETSc Viewer is used for writing PETSc Mat object in binary format */
PetscViewer view;
/*Just record time required for conversion */
PetscLogDouble t1,t2,elapsed_time;
/* MatrixMarket struct */
MM_typecode matcode;
/*Initialise PETSc lib */
PetscInitialize(&argc,&args,(char *)0,PETSC_NULL);
/* Just record time */
//ierr = PetscGetTime(&t1); CHKERRQ(ierr);
/*Get name of matrix market file from command line options and Open file*/
ierr = PetscOptionsGetString(PETSC_NULL,"-fin",inputFile,127,PETSC_NULL); CHKERRQ(ierr);
ierr = PetscFOpen(PETSC_COMM_SELF,inputFile,"r",&file); CHKERRQ(ierr);
if (mm_read_banner(file, &matcode)) {
PetscPrintf(PETSC_COMM_SELF, "Could not read Matrix Market banner.\n");
exit(1);
}
/********************* MM_typecode query fucntions ***************************/
/* #define mm_is_matrix(typecode) ((typecode)[0]=='M') */
/* #define mm_is_sparse(typecode) ((typecode)[1]=='C') */
/* #define mm_is_coordinate(typecode)((typecode)[1]=='C') */
/* #define mm_is_dense(typecode) ((typecode)[1]=='A') */
/* #define mm_is_array(typecode) ((typecode)[1]=='A') */
/* #define mm_is_complex(typecode) ((typecode)[2]=='C') */
/* #define mm_is_real(typecode) ((typecode)[2]=='R') */
/* #define mm_is_pattern(typecode) ((typecode)[2]=='P') */
/* #define mm_is_integer(typecode) ((typecode)[2]=='I') */
/* #define mm_is_symmetric(typecode)((typecode)[3]=='S') */
/* #define mm_is_general(typecode) ((typecode)[3]=='G') */
/* #define mm_is_skew(typecode) ((typecode)[3]=='K') */
/* #define mm_is_hermitian(typecode)((typecode)[3]=='H') */
/* int mm_is_valid(MM_typecode matcode); */
/* Do not convert pattern matrices */
if (mm_is_pattern(matcode)) {
ierr = PetscPrintf(PETSC_COMM_SELF, "%s: Pattern matrix -- skipping.\n", inputFile);
exit(0);
}
/* find out size of sparse matrix .... */
/*Reads size of sparse matrix from matrix market file */
int ret_code;
if ((ret_code = mm_read_mtx_crd_size(file, &m, &n, &nnz)) !=0)
exit(1);
ierr = PetscPrintf(PETSC_COMM_SELF, "%s: ROWS = %d, COLUMNS = %d, NO OF NON-ZEROS = %d\n",inputFile,m,n,nnz);
/* Only consider square matrices */
if (m != n) {
ierr = PetscPrintf(PETSC_COMM_SELF, "%s: Nonsquare matrix -- skipping.\n", inputFile);
exit(0);
}
ierr = MatCreate(PETSC_COMM_WORLD,&pMat);CHKERRQ(ierr);
ierr = MatSetFromOptions(pMat);CHKERRQ(ierr);
//ierr = MatSetOption(pMat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); CHKERRQ(ierr);
if (mm_is_symmetric(matcode)) {
ierr = MatSetOption(pMat,MAT_SYMMETRIC,PETSC_TRUE); CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:kanika901,项目名称:lighthouse,代码行数:101,代码来源:MatMarket_to_PETScBin.c
示例10: main
PetscInt main(PetscInt argc,char **args)
{
Mat A,As;
PetscBool flg,disp_mat=PETSC_FALSE;
PetscErrorCode ierr;
PetscMPIInt size,rank;
PetscInt i,j;
PetscScalar v,sigma2;
PetscRandom rctx;
PetscReal h2,sigma1=100.0;
PetscInt dim,Ii,J,n = 3,use_random,rstart,rend;
KSP ksp;
PC pc;
Mat F;
PetscInt nneg, nzero, npos;
PetscInitialize(&argc,&args,(char *)0,help);
#if !defined(PETSC_USE_COMPLEX)
SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
#endif
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = PetscOptionsHasName(PETSC_NULL, "-display_mat", &disp_mat);CHKERRQ(ierr);
ierr = PetscOptionsGetReal(PETSC_NULL,"-sigma1",&sigma1,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
dim = n*n;
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
ierr = MatSetType(A,MATAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = PetscOptionsHasName(PETSC_NULL,"-norandom",&flg);CHKERRQ(ierr);
if (flg) use_random = 0;
else use_random = 1;
if (use_random) {
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
ierr = PetscRandomSetInterval(rctx,0.0,PETSC_i);CHKERRQ(ierr);
ierr = PetscRandomGetValue(rctx,&sigma2);CHKERRQ(ierr); /* RealPart(sigma2) == 0.0 */
} else {
sigma2 = 10.0*PETSC_i;
}
h2 = 1.0/((n+1)*(n+1));
ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr);
for (Ii=rstart; Ii<rend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {
J = Ii-n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (i<n-1) {
J = Ii+n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j>0) {
J = Ii-1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j<n-1) {
J = Ii+1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
v = 4.0 - sigma1*h2;
ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Check whether A is symmetric */
ierr = PetscOptionsHasName(PETSC_NULL, "-check_symmetric", &flg);CHKERRQ(ierr);
if (flg) {
Mat Trans;
ierr = MatTranspose(A,MAT_INITIAL_MATRIX, &Trans);
ierr = MatEqual(A, Trans, &flg);
if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"A is not symmetric");
ierr = MatDestroy(&Trans);CHKERRQ(ierr);
}
ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
/* make A complex Hermitian */
Ii = 0; J = dim-1;
if (Ii >= rstart && Ii < rend){
v = sigma2*h2; /* RealPart(v) = 0.0 */
ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
v = -sigma2*h2;
ierr = MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
Ii = dim-2; J = dim-1;
if (Ii >= rstart && Ii < rend){
v = sigma2*h2; /* RealPart(v) = 0.0 */
ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);
v = -sigma2*h2;
ierr = MatSetValues(A,1,&J,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* Check whether A is Hermitian */
ierr = PetscOptionsHasName(PETSC_NULL, "-check_Hermitian", &flg);CHKERRQ(ierr);
if (flg) {
Mat Hermit;
if (disp_mat){
if (!rank) printf(" A:\n");
//.........这里部分代码省略.........
开发者ID:erdc-cm,项目名称:petsc-dev,代码行数:101,代码来源:ex36.c
示例11: main
int main(int argc,char **args)
{
Mat A,B,F;
PetscErrorCode ierr;
KSP ksp;
PC pc;
PetscInt N, n=10, m, Istart, Iend, II, J, i,j;
PetscInt nneg, nzero, npos;
PetscScalar v,sigma;
PetscBool flag,loadA=PETSC_FALSE,loadB=PETSC_FALSE;
char file[2][PETSC_MAX_PATH_LEN];
PetscViewer viewer;
PetscMPIInt rank;
PetscInitialize(&argc,&args,(char *)0,help);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrices that define the eigensystem, Ax=kBx
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = PetscOptionsGetString(PETSC_NULL,"-fA",file[0],PETSC_MAX_PATH_LEN,&loadA);CHKERRQ(ierr);
if (loadA) {
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[0],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetType(A,MATSBAIJ);CHKERRQ(ierr);
ierr = MatLoad(A,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = PetscOptionsGetString(PETSC_NULL,"-fB",file[1],PETSC_MAX_PATH_LEN,&loadB);CHKERRQ(ierr);
if (loadB){
/* load B to get A = A + sigma*B */
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[1],FILE_MODE_READ,&viewer);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
ierr = MatSetType(B,MATSBAIJ);CHKERRQ(ierr);
ierr = MatLoad(B,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
}
}
if (!loadA) { /* Matrix A is copied from slepc-3.0.0-p6/src/examples/ex13.c. */
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);CHKERRQ(ierr);
if( flag==PETSC_FALSE ) m=n;
N = n*m;
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr);
ierr = MatSetType(A,MATSBAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = MatSetOption(A,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
for( II=Istart; II<Iend; II++ ) {
v = -1.0; i = II/n; j = II-i*n;
if(i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr); }
if(i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr); }
if(j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr); }
if(j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr); }
v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
/* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
if (!loadB) {
ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
ierr = MatSetSizes(B,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
ierr = MatSetType(B,MATSBAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(B);CHKERRQ(ierr);
ierr = MatSetUp(B);CHKERRQ(ierr);
ierr = MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
for( II=Istart; II<Iend; II++ ) {
/* v=4.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);CHKERRQ(ierr); */
v=1.0; MatSetValues(B,1,&II,1,&II,&v,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
/* ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
/* Set a shift: A = A - sigma*B */
ierr = PetscOptionsGetScalar(PETSC_NULL,"-sigma",&sigma,&flag);CHKERRQ(ierr);
if (flag){
sigma = -1.0 * sigma;
ierr = MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); /* A <- A - sigma*B */
/* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
}
/* Test MatGetInertia() */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPPREONLY);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCCHOLESKY);CHKERRQ(ierr);
ierr = PCSetFromOptions(pc);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:Kun-Qu,项目名称:petsc,代码行数:101,代码来源:ex33.c
-
六六分期app的软件客服如何联系?不知道吗?加qq群【895510560】即可!标题:六六分期
阅读:19247|2023-10-27
-
今天小编告诉大家如何处理win10系统火狐flash插件总是崩溃的问题,可能很多用户都不知
阅读:10004|2022-11-06
-
今天小编告诉大家如何对win10系统删除桌面回收站图标进行设置,可能很多用户都不知道
阅读:8334|2022-11-06
-
今天小编告诉大家如何对win10系统电脑设置节能降温的设置方法,想必大家都遇到过需要
阅读:8703|2022-11-06
-
我们在使用xp系统的过程中,经常需要对xp系统无线网络安装向导设置进行设置,可能很多
阅读:8649|2022-11-06
-
今天小编告诉大家如何处理win7系统玩cf老是与主机连接不稳定的问题,可能很多用户都不
阅读:9675|2022-11-06
-
电脑对日常生活的重要性小编就不多说了,可是一旦碰到win7系统设置cf烟雾头的问题,很
阅读:8635|2022-11-06
-
我们在日常使用电脑的时候,有的小伙伴们可能在打开应用的时候会遇见提示应用程序无法
阅读:8008|2022-11-06
-
今天小编告诉大家如何对win7系统打开vcf文件进行设置,可能很多用户都不知道怎么对win
阅读:8671|2022-11-06
-
今天小编告诉大家如何对win10系统s4开启USB调试模式进行设置,可能很多用户都不知道怎
阅读:7542|2022-11-06
|
请发表评论