本文整理汇总了C++中MatSetSizes函数的典型用法代码示例。如果您正苦于以下问题：C++ MatSetSizes函数的具体用法？C++ MatSetSizes怎么用？C++ MatSetSizes使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了MatSetSizes函数的20个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于我们的系统推荐出更棒的C++代码示例。
示例1: main
int main(int argc,char **args)
{
PetscMPIInt size;
PetscErrorCode ierr;
Vec x,y,b,s1,s2;
Mat A; /* linear system matrix */
Mat sA,sB,sC; /* symmetric part of the matrices */
PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc;
PetscReal norm1,norm2,rnorm,tol=PETSC_SMALL;
PetscScalar neg_one = -1.0,four=4.0,value[3];
IS perm, iscol;
PetscRandom rdm;
PetscBool doIcc=PETSC_TRUE,equal;
MatInfo minfo1,minfo2;
MatFactorInfo factinfo;
MatType type;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
ierr = PetscOptionsGetInt(NULL,"-bs",&bs,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr);
n = mbs*bs;
ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSetType(A,MATSEQBAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSeqBAIJSetPreallocation(A,bs,nz,NULL);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_SELF,&sA);CHKERRQ(ierr);
ierr = MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSetType(sA,MATSEQSBAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(sA);CHKERRQ(ierr);
ierr = MatGetType(sA,&type);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc);CHKERRQ(ierr);
ierr = MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL);CHKERRQ(ierr);
ierr = MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr);
/* Test MatGetOwnershipRange() */
ierr = MatGetOwnershipRange(A,&Ii,&J);CHKERRQ(ierr);
ierr = MatGetOwnershipRange(sA,&i,&j);CHKERRQ(ierr);
if (i-Ii || j-J) {
ierr = PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n");CHKERRQ(ierr);
}
/* Assemble matrix */
if (bs == 1) {
ierr = PetscOptionsGetInt(NULL,"-test_problem",&prob,NULL);CHKERRQ(ierr);
if (prob == 1) { /* tridiagonal matrix */
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=1; i<n-1; i++) {
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1;
value[0]= 0.1; value[1]=-1; value[2]=2;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
i = 0;
col[0] = n-1; col[1] = 1; col[2] = 0;
value[0] = 0.1; value[1] = -1.0; value[2] = 2;
ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
} else if (prob ==2) { /* matrix for the five point stencil */
n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001);
if (n1*n1 - n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive interger!");
for (i=0; i<n1; i++) {
for (j=0; j<n1; j++) {
Ii = j + n1*i;
if (i>0) {
J = Ii - n1;
ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
}
if (i<n1-1) {
J = Ii + n1;
ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
}
if (j>0) {
J = Ii - 1;
ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
}
if (j<n1-1) {
J = Ii + 1;
ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr);
}
}
//.........这里部分代码省略.........
开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:ex74.c
示例2: main
int main(int argc,char **argv)
{
PetscErrorCode ierr;
Mat seqmat,mpimat;
PetscMPIInt rank;
PetscScalar value[3],*vals;
PetscInt i,col[3],n=5,bs=1;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-bs",&bs,NULL);CHKERRQ(ierr);
/* Create seqaij matrices of size (n+rank) by n */
ierr = MatCreate(PETSC_COMM_SELF,&seqmat);CHKERRQ(ierr);
ierr = MatSetSizes(seqmat,(n+rank)*bs,PETSC_DECIDE,PETSC_DECIDE,n*bs);CHKERRQ(ierr);
ierr = MatSetFromOptions(seqmat);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(seqmat,3*bs,NULL);CHKERRQ(ierr);
ierr = MatSeqBAIJSetPreallocation(seqmat,bs,3,NULL);CHKERRQ(ierr);
if (bs == 1) {
value[0] = -1.0; value[1] = 2.0; value[2] = -1.0;
for (i=1; i<n-1; i++) {
col[0] = i-1; col[1] = i; col[2] = i+1;
ierr = MatSetValues(seqmat,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr);
}
i = n - 1; col[0] = n - 2; col[1] = n - 1;
ierr = MatSetValues(seqmat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
i = 0; col[0] = 0; col[1] = 1; value[0] = 2.0; value[1] = -1.0;
ierr = MatSetValues(seqmat,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr);
} else {
PetscInt *rows,*cols,j;
ierr = PetscMalloc3(bs*bs,&vals,bs,&rows,bs,&cols);CHKERRQ(ierr);
/* diagonal blocks */
for (i=0; i<bs*bs; i++) vals[i] = 2.0;
for (i=0; i<n*bs; i+=bs) {
for (j=0; j<bs; j++) {rows[j] = i+j; cols[j] = i+j;}
ierr = MatSetValues(seqmat,bs,rows,bs,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
}
/* off-diagonal blocks */
for (i=0; i<bs*bs; i++) vals[i] = -1.0;
for (i=0; i<(n-1)*bs; i+=bs) {
for (j=0; j<bs; j++) {rows[j] = i+j; cols[j] = i+bs+j;}
ierr = MatSetValues(seqmat,bs,rows,bs,cols,vals,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = PetscFree3(vals,rows,cols);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(seqmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(seqmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (rank == 1) {
printf("[%d] seqmat:\n",rank);
ierr = MatView(seqmat,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
}
ierr = MatCreateMPIMatConcatenateSeqMat(PETSC_COMM_WORLD,seqmat,PETSC_DECIDE,MAT_INITIAL_MATRIX,&mpimat);CHKERRQ(ierr);
ierr = MatCreateMPIMatConcatenateSeqMat(PETSC_COMM_WORLD,seqmat,PETSC_DECIDE,MAT_REUSE_MATRIX,&mpimat);CHKERRQ(ierr);
ierr = MatView(mpimat,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = MatDestroy(&seqmat);CHKERRQ(ierr);
ierr = MatDestroy(&mpimat);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:fengyuqi，项目名称:petsc，代码行数:64，代码来源:ex86.c
示例3: main
int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* KSP context */
PetscErrorCode ierr;
PetscInt n = 10,its, dim,p = 1,use_random;
PetscScalar none = -1.0,pfive = 0.5;
PetscReal norm;
PetscRandom rctx;
TestType type;
PetscBool flg;
PetscInitialize(&argc,&args,(char *)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);
switch (p) {
case 1: type = TEST_1; dim = n; break;
case 2: type = TEST_2; dim = n; break;
case 3: type = TEST_3; dim = n; break;
case 4: type = HELMHOLTZ_1; dim = n*n; break;
case 5: type = HELMHOLTZ_2; dim = n*n; break;
default: type = TEST_1; dim = n;
}
/* Create vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,dim);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
use_random = 1;
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-norandom",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
use_random = 0;
ierr = VecSet(u,pfive);CHKERRQ(ierr);
} else {
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
}
/* Create and assemble matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = FormTestMatrix(A,n,type);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-printout",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/* Create KSP context; set operators and options; solve linear system */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Check error */
ierr = VecAXPY(x,none,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G,Iterations %D\n",norm,its);CHKERRQ(ierr);
/* Free work space */
ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);}
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:erdc-cm，项目名称:petsc-dev，代码行数:80，代码来源:ex17.c
示例4: main
int main(int argc,char **args)
{
Mat C;
PetscScalar v,none = -1.0;
PetscInt i,j,Ii,J,Istart,Iend,N,m = 4,n = 4,its,k;
PetscErrorCode ierr;
PetscMPIInt size,rank;
PetscReal err_norm,res_norm,err_tol=1.e-7,res_tol=1.e-6;
Vec x,b,u,u_tmp;
PetscRandom r;
PC pc;
KSP ksp;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
N = m*n;
/* Generate 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);
ierr = MatGetOwnershipRange(C,&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; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
v = 4.0; ierr = MatSetValues(C,1,&Ii,1,&Ii,&v,ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/* a shift can make C indefinite. Preconditioners LU, ILU (for BAIJ format) and ICC may fail */
/* ierr = MatShift(C,alpha);CHKERRQ(ierr); */
/* ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
/* Setup and solve for system */
/* Create vectors. */
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u_tmp);CHKERRQ(ierr);
/* Set exact solution u; then compute right-hand-side vector b. */
ierr = PetscRandomCreate(PETSC_COMM_SELF,&r);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(r);CHKERRQ(ierr);
ierr = VecSetRandom(u,r);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
ierr = MatMult(C,u,b);CHKERRQ(ierr);
for (k=0; k<3; k++) {
if (k == 0) { /* CG */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,C,C);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n CG: \n");CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
} else if (k == 1) { /* MINRES */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,C,C);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n MINRES: \n");CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPMINRES);CHKERRQ(ierr);
} else { /* SYMMLQ */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,C,C);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n SYMMLQ: \n");CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPSYMMLQ);CHKERRQ(ierr);
}
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
/* ierr = PCSetType(pc,PCICC);CHKERRQ(ierr); */
ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* Solve linear system; */
ierr = KSPSetUp(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
/* Check error */
ierr = VecCopy(u,u_tmp);CHKERRQ(ierr);
ierr = VecAXPY(u_tmp,none,x);CHKERRQ(ierr);
ierr = VecNorm(u_tmp,NORM_2,&err_norm);CHKERRQ(ierr);
ierr = MatMult(C,x,u_tmp);CHKERRQ(ierr);
ierr = VecAXPY(u_tmp,none,b);CHKERRQ(ierr);
ierr = VecNorm(u_tmp,NORM_2,&res_norm);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:ex24.c
示例5: main
Modified from the code contributed by Yaning Liu @lbl.gov \n\n";
/*
Example:
mpiexec -n <np> ./ex103
mpiexec -n <np> ./ex103 -mat_type elemental -mat_view
mpiexec -n <np> ./ex103 -mat_type aij
*/
#include <petscmat.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char** argv)
{
Mat A,A_elemental;
PetscInt i,j,M=10,N=5,nrows,ncols;
PetscErrorCode ierr;
PetscMPIInt rank,size;
IS isrows,iscols;
const PetscInt *rows,*cols;
PetscScalar *v;
MatType type;
PetscBool isDense,isAIJ,flg;
ierr = PetscInitialize(&argc, &argv, (char*)0, help);if (ierr) return ierr;
#if !defined(PETSC_HAVE_ELEMENTAL)
SETERRQ(PETSC_COMM_WORLD,1,"This example requires ELEMENTAL");
#endif
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
/* Creat a matrix */
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);
/* 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++) {
if (size == 1) {
v[i*ncols+j] = (PetscScalar)(i+j);
} else {
v[i*ncols+j] = (PetscScalar)rank+j*0.1;
}
}
}
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 = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%D] local nrows %D, ncols %D\n",rank,nrows,ncols);CHKERRQ(ierr);
//ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
/* Test MatSetValues() by converting A to A_elemental */
ierr = MatGetType(A,&type);CHKERRQ(ierr);
if (size == 1) {
ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQDENSE,&isDense);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&isAIJ);CHKERRQ(ierr);
} else {
ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIDENSE,&isDense);CHKERRQ(ierr);
ierr = PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&isAIJ);CHKERRQ(ierr);
}
if (isDense || isAIJ) {
Mat Aexplicit;
ierr = MatConvert(A, MATELEMENTAL, MAT_INITIAL_MATRIX, &A_elemental);CHKERRQ(ierr);
ierr = MatComputeExplicitOperator(A_elemental,&Aexplicit);CHKERRQ(ierr);
ierr = MatMultEqual(Aexplicit,A_elemental,5,&flg);CHKERRQ(ierr);
if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Aexplicit != A_elemental.");
ierr = MatDestroy(&Aexplicit);CHKERRQ(ierr);
/* Test MAT_REUSE_MATRIX which is only supported for inplace conversion */
ierr = MatConvert(A, MATELEMENTAL, MAT_INPLACE_MATRIX, &A);CHKERRQ(ierr);
ierr = MatMultEqual(A_elemental,A,5,&flg);CHKERRQ(ierr);
if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A_elemental != A.");
ierr = MatDestroy(&A_elemental);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 = PetscFree(v);CHKERRQ(ierr);
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
开发者ID:tom-klotz，项目名称:petsc，代码行数:97，代码来源:ex103.c
示例6: main
int main(int argc,char **args)
{
Mat C;
int i,m = 5,rank,size,N,start,end,M;
int ierr,idx[4];
PetscBool flg;
PetscScalar Ke[16];
PetscReal h;
Vec u,b;
KSP ksp;
MatNullSpace nullsp;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = PetscOptionsGetInt(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);
start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
end = start + M/size + ((M%size) > rank);
/* Assemble matrix */
ierr = FormElementStiffness(h*h,Ke); /* element stiffness for Laplacian */
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 right-hand-side and solution vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)u,"Approx. Solution");CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)b,"Right hand side");CHKERRQ(ierr);
ierr = VecSet(u,1.0);CHKERRQ(ierr);
ierr = MatMult(C,u,b);CHKERRQ(ierr);
ierr = VecSet(u,0.0);CHKERRQ(ierr);
/* Solve linear system */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(NULL,"-fixnullspace",&flg,NULL);CHKERRQ(ierr);
if (flg) {
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,NULL,&nullsp);CHKERRQ(ierr);
ierr = KSPSetNullSpace(ksp,nullsp);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&&nullsp);CHKERRQ(ierr);
}
ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);
/* Free work space */
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
开发者ID:feelpp，项目名称:debian-petsc，代码行数:75，代码来源:ex20.c
示例7: g2_correlation
void g2_correlation(PetscScalar ***g2_values,Vec dm0,PetscInt n_tau,PetscReal tau_max,PetscInt n_st,PetscReal st_max,PetscInt number_of_ops,...){
TSCtx tsctx;
PetscReal st_dt,previous_start_time,this_start_time,dt;
PetscReal tau_t_max,dt_tau;
PetscInt i,j,dim,steps_max;
Mat A_star_A,tmp_mat;
Vec init_dm;
va_list ap;
/*Explicitly construct our jump matrix by adding up all of the operators
* \rho = A \rho A^\dag
* Vectorized:
* \rho = (A* \cross A) \rho
*/
dim = total_levels*total_levels; //Assumes Lindblad
MatCreate(PETSC_COMM_WORLD,&tmp_mat);
MatSetType(tmp_mat,MATMPIAIJ);
MatSetSizes(tmp_mat,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
MatSetFromOptions(tmp_mat);
MatMPIAIJSetPreallocation(tmp_mat,4,NULL,4,NULL);
va_start(ap,number_of_ops);
//Get A* \cross I
vadd_ops_to_mat(tmp_mat,1,number_of_ops,ap);
va_end(ap);
MatCreate(PETSC_COMM_WORLD,&tsctx.I_cross_A);
MatSetType(tsctx.I_cross_A,MATMPIAIJ);
MatSetSizes(tsctx.I_cross_A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);
MatSetFromOptions(tsctx.I_cross_A);
MatMPIAIJSetPreallocation(tsctx.I_cross_A,4,NULL,4,NULL);
va_start(ap,number_of_ops);
//Get I_cross_A
vadd_ops_to_mat(tsctx.I_cross_A,-1,number_of_ops,ap);
va_end(ap);
//Get (A* \cross I) (I \cross A)
MatMatMult(tmp_mat,tsctx.I_cross_A,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&A_star_A);
MatDestroy(&tmp_mat);
VecDuplicate(dm0,&(tsctx.tmp_dm));
VecDuplicate(dm0,&(tsctx.tmp_dm2));
VecDuplicate(dm0,&init_dm);
//Arbitrary, set this better
steps_max = 51000;
dt = 0.025;
//Allocate memory for g2
(*g2_values) = (PetscScalar **)malloc((n_st+1)*sizeof(PetscScalar *));
for (i=0;i<n_st+1;i++){
(*g2_values)[i] = (PetscScalar *)malloc((n_tau+1)*sizeof(PetscScalar));
for (j=0;j<n_tau+1;j++){
(*g2_values)[i][j] = 0.0;
}
}
tsctx.g2_values = (*g2_values);
set_ts_monitor_ctx(_g2_ts_monitor,&tsctx);
st_dt = st_max/n_st;
previous_start_time = 0;
tsctx.i_st = 0;
tsctx.i_st = tsctx.i_st + 1; //Why?
for (this_start_time=st_dt;this_start_time<=st_max;this_start_time+=st_dt){
//Go from previous start time to this_start_time
tsctx.tau_evolve = 0;
dt = (this_start_time - previous_start_time)/500; //500 is arbitrary, should be picked better
time_step(dm0,previous_start_time,this_start_time,dt,steps_max);
//Timestep through taus
tau_t_max = this_start_time + tau_max;
dt_tau = (tau_t_max - this_start_time)/n_tau;
/*
* Force an 'emission' to get A \rho A^\dag terms
* We already have A* \cross A - we just do the multiplication
*/
tsctx.tau_evolve = 1;
//Copy the timestepped dm into our init_dm for tau sweep
VecCopy(dm0,init_dm);
MatMult(A_star_A,dm0,init_dm); //init_dm = A * dm0
tsctx.i_tau = 0;
time_step(init_dm,this_start_time,tau_t_max,dt_tau,steps_max);
previous_start_time = this_start_time;
tsctx.i_st = tsctx.i_st + 1;
}
return;
}
开发者ID:0tt3r，项目名称:QuaC，代码行数:92，代码来源:solver.c
示例8: PCGAMGProlongator_Classical_Direct
//.........这里部分代码省略.........
lsparse[i] = 0;
gsparse[i] = 0;
if (lcid[i] >= 0) {
lsparse[i] = 1;
gsparse[i] = 0;
} else {
ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
for (j = 0;j < ncols;j++) {
col = rcol[j];
if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
lsparse[i] += 1;
}
}
ierr = MatRestoreRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
/* off */
if (gA) {
ierr = MatGetRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
for (j = 0; j < ncols; j++) {
col = rcol[j];
if (gcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
gsparse[i] += 1;
}
}
ierr = MatRestoreRow(gA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
}
}
}
/* preallocate and create the prolongator */
ierr = MatCreate(comm,P); CHKERRQ(ierr);
ierr = MatGetType(G,&mtype);CHKERRQ(ierr);
ierr = MatSetType(*P,mtype);CHKERRQ(ierr);
ierr = MatSetSizes(*P,fn,cn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(*P,0,lsparse,0,gsparse);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(*P,0,lsparse);CHKERRQ(ierr);
/* loop over local fine nodes -- get the diagonal, the sum of positive and negative strong and weak weights, and set up the row */
for (i = 0;i < fn;i++) {
/* determine on or off */
row_f = i + fs;
row_c = lcid[i];
if (row_c >= 0) {
pij = 1.;
ierr = MatSetValues(*P,1,&row_f,1,&row_c,&pij,INSERT_VALUES);CHKERRQ(ierr);
} else {
g_pos = 0.;
g_neg = 0.;
a_pos = 0.;
a_neg = 0.;
diag = 0.;
/* local connections */
ierr = MatGetRow(lA,i,&ncols,&rcol,&rval);CHKERRQ(ierr);
for (j = 0; j < ncols; j++) {
col = rcol[j];
if (lcid[col] >= 0 && (PetscRealPart(rval[j]) > gamg->threshold*Amax_pos[i] || PetscRealPart(-rval[j]) > gamg->threshold*Amax_neg[i])) {
if (PetscRealPart(rval[j]) > 0.) {
g_pos += rval[j];
} else {
g_neg += rval[j];
}
}
if (col != i) {
if (PetscRealPart(rval[j]) > 0.) {
a_pos += rval[j];
开发者ID:ZJLi2013，项目名称:petsc，代码行数:67，代码来源:classical.c
示例9: PCGAMGTruncateProlongator_Private
PetscErrorCode PCGAMGTruncateProlongator_Private(PC pc,Mat *P)
{
PetscInt j,i,ps,pf,pn,pcs,pcf,pcn,idx,cmax;
PetscErrorCode ierr;
const PetscScalar *pval;
const PetscInt *pcol;
PetscScalar *pnval;
PetscInt *pncol;
PetscInt ncols;
Mat Pnew;
PetscInt *lsparse,*gsparse;
PetscReal pmax_pos,pmax_neg,ptot_pos,ptot_neg,pthresh_pos,pthresh_neg;
PC_MG *mg = (PC_MG*)pc->data;
PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
PC_GAMG_Classical *cls = (PC_GAMG_Classical*)pc_gamg->subctx;
PetscFunctionBegin;
/* trim and rescale with reallocation */
ierr = MatGetOwnershipRange(*P,&ps,&pf);CHKERRQ(ierr);
ierr = MatGetOwnershipRangeColumn(*P,&pcs,&pcf);CHKERRQ(ierr);
pn = pf-ps;
pcn = pcf-pcs;
ierr = PetscMalloc(sizeof(PetscInt)*pn,&lsparse);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscInt)*pn,&gsparse);CHKERRQ(ierr);
/* allocate */
cmax = 0;
for (i=ps;i<pf;i++) {
lsparse[i-ps] = 0;
gsparse[i-ps] = 0;
ierr = MatGetRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
if (ncols > cmax) {
cmax = ncols;
}
pmax_pos = 0.;
pmax_neg = 0.;
for (j=0;j<ncols;j++) {
if (PetscRealPart(pval[j]) > pmax_pos) {
pmax_pos = PetscRealPart(pval[j]);
} else if (PetscRealPart(pval[j]) < pmax_neg) {
pmax_neg = PetscRealPart(pval[j]);
}
}
for (j=0;j<ncols;j++) {
if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold || PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
if (pcol[j] >= pcs && pcol[j] < pcf) {
lsparse[i-ps]++;
} else {
gsparse[i-ps]++;
}
}
}
ierr = MatRestoreRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
}
ierr = PetscMalloc(sizeof(PetscScalar)*cmax,&pnval);CHKERRQ(ierr);
ierr = PetscMalloc(sizeof(PetscInt)*cmax,&pncol);CHKERRQ(ierr);
ierr = MatCreate(PetscObjectComm((PetscObject)*P),&Pnew);CHKERRQ(ierr);
ierr = MatSetType(Pnew, MATAIJ);CHKERRQ(ierr);
ierr = MatSetSizes(Pnew,pn,pcn,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(Pnew,0,lsparse);CHKERRQ(ierr);
ierr = MatMPIAIJSetPreallocation(Pnew,0,lsparse,0,gsparse);CHKERRQ(ierr);
for (i=ps;i<pf;i++) {
ierr = MatGetRow(*P,i,&ncols,&pcol,&pval);CHKERRQ(ierr);
pmax_pos = 0.;
pmax_neg = 0.;
for (j=0;j<ncols;j++) {
if (PetscRealPart(pval[j]) > pmax_pos) {
pmax_pos = PetscRealPart(pval[j]);
} else if (PetscRealPart(pval[j]) < pmax_neg) {
pmax_neg = PetscRealPart(pval[j]);
}
}
pthresh_pos = 0.;
pthresh_neg = 0.;
ptot_pos = 0.;
ptot_neg = 0.;
for (j=0;j<ncols;j++) {
if (PetscRealPart(pval[j]) >= cls->interp_threshold*pmax_pos) {
pthresh_pos += PetscRealPart(pval[j]);
} else if (PetscRealPart(pval[j]) <= cls->interp_threshold*pmax_neg) {
pthresh_neg += PetscRealPart(pval[j]);
}
if (PetscRealPart(pval[j]) > 0.) {
ptot_pos += PetscRealPart(pval[j]);
} else {
ptot_neg += PetscRealPart(pval[j]);
}
}
if (PetscAbsReal(pthresh_pos) > 0.) ptot_pos /= pthresh_pos;
if (PetscAbsReal(pthresh_neg) > 0.) ptot_neg /= pthresh_neg;
idx=0;
for (j=0;j<ncols;j++) {
if (PetscRealPart(pval[j]) >= pmax_pos*cls->interp_threshold) {
pnval[idx] = ptot_pos*pval[j];
pncol[idx] = pcol[j];
idx++;
} else if (PetscRealPart(pval[j]) <= pmax_neg*cls->interp_threshold) {
pnval[idx] = ptot_neg*pval[j];
//.........这里部分代码省略.........
开发者ID:ZJLi2013，项目名称:petsc，代码行数:101，代码来源:classical.c
示例10: main
int main(int argc,char **argv)
{
SNES snes; /* SNES context */
Vec x,r,F,U; /* vectors */
Mat J; /* Jacobian matrix */
MonitorCtx monP; /* monitoring context */
PetscErrorCode ierr;
PetscInt its,n = 5,i,maxit,maxf;
PetscMPIInt size;
PetscScalar h,xp,v,none = -1.0;
PetscReal abstol,rtol,stol,norm;
PetscInitialize(&argc,&argv,(char*)0,help);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
h = 1.0/(n-1);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create nonlinear solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create vector data structures; set function evaluation routine
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Note that we form 1 vector from scratch and then duplicate as needed.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
ierr = VecDuplicate(x,&F);CHKERRQ(ierr);
ierr = VecDuplicate(x,&U);CHKERRQ(ierr);
/*
Set function evaluation routine and vector
*/
ierr = SNESSetFunction(snes,r,FormFunction,(void*)F);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create matrix data structure; set Jacobian evaluation routine
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(J);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(J,3,NULL);CHKERRQ(ierr);
/*
Set Jacobian matrix data structure and default Jacobian evaluation
routine. User can override with:
-snes_fd : default finite differencing approximation of Jacobian
-snes_mf : matrix-free Newton-Krylov method with no preconditioning
(unless user explicitly sets preconditioner)
-snes_mf_operator : form preconditioning matrix as set by the user,
but use matrix-free approx for Jacobian-vector
products within Newton-Krylov method
*/
ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Customize nonlinear solver; set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Set an optional user-defined monitoring routine
*/
ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);CHKERRQ(ierr);
ierr = SNESMonitorSet(snes,Monitor,&monP,0);CHKERRQ(ierr);
/*
Set names for some vectors to facilitate monitoring (optional)
*/
ierr = PetscObjectSetName((PetscObject)x,"Approximate Solution");CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject)U,"Exact Solution");CHKERRQ(ierr);
/*
Set SNES/KSP/KSP/PC runtime options, e.g.,
-snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
*/
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
/*
Print parameters used for convergence testing (optional) ... just
to demonstrate this routine; this information is also printed with
the option -snes_view
*/
ierr = SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Initialize application:
Store right-hand-side of PDE and exact solution
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//.........这里部分代码省略.........
开发者ID:00liujj，项目名称:petsc，代码行数:101，代码来源:ex2.c
示例11: main
int main(int argc,char **argv)
{
TS ts; /* nonlinear solver */
Vec x; /* solution, residual vectors */
Mat A; /* Jacobian matrix */
PetscInt steps;
PetscReal ftime = 0.5;
PetscBool monitor = PETSC_FALSE;
PetscScalar *x_ptr;
PetscMPIInt size;
struct _n_User user;
PetscErrorCode ierr;
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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!");
ierr = RegisterMyARK2();CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set runtime options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
user.next_output = 0.0;
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);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create timestepping solver context
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);
ierr = TSSetIJacobian(ts,A,A,IJacobian,&user);CHKERRQ(ierr);
ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime);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] = -2.355301397608119909925287735864250951918;
ierr = VecRestoreArray(x,&x_ptr);CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,0.0,.001);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,"steps %D, ftime %g\n",steps,(double)ftime);CHKERRQ(ierr);
ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = PetscFinalize();
PetscFunctionReturn(0);
}
开发者ID:wgapl，项目名称:petsc，代码行数:84，代码来源:ex19.c
示例12: solvers
EXTERN_C_END
/*MC
MATSOLVERUMFPACK = "umfpack" - A matrix type providing direct solvers (LU) for sequential matrices
via the external package UMFPACK.
./configure --download-umfpack to install PETSc to use UMFPACK
Consult UMFPACK documentation for more information about the Control parameters
which correspond to the options database keys below.
Options Database Keys:
+ -mat_umfpack_prl - UMFPACK print level: Control[UMFPACK_PRL]
. -mat_umfpack_strategy <AUTO> - (choose one of) AUTO UNSYMMETRIC SYMMETRIC 2BY2
. -mat_umfpack_dense_col <alpha_c> - UMFPACK dense column threshold: Control[UMFPACK_DENSE_COL]
. -mat_umfpack_dense_row <0.2> - Control[UMFPACK_DENSE_ROW]
. -mat_umfpack_amd_dense <10> - Control[UMFPACK_AMD_DENSE]
. -mat_umfpack_block_size <bs> - UMFPACK block size for BLAS-Level 3 calls: Control[UMFPACK_BLOCK_SIZE]
. -mat_umfpack_2by2_tolerance <0.01> - Control[UMFPACK_2BY2_TOLERANCE]
. -mat_umfpack_fixq <0> - Control[UMFPACK_FIXQ]
. -mat_umfpack_aggressive <1> - Control[UMFPACK_AGGRESSIVE]
. -mat_umfpack_pivot_tolerance <delta> - UMFPACK partial pivot tolerance: Control[UMFPACK_PIVOT_TOLERANCE]
. -mat_umfpack_sym_pivot_tolerance <0.001> - Control[UMFPACK_SYM_PIVOT_TOLERANCE]
. -mat_umfpack_scale <NONE> - (choose one of) NONE SUM MAX
. -mat_umfpack_alloc_init <delta> - UMFPACK factorized matrix allocation modifier: Control[UMFPACK_ALLOC_INIT]
. -mat_umfpack_droptol <0> - Control[UMFPACK_DROPTOL]
- -mat_umfpack_irstep <maxit> - UMFPACK maximum number of iterative refinement steps: Control[UMFPACK_IRSTEP]
Level: beginner
.seealso: PCLU, MATSOLVERSUPERLU, MATSOLVERMUMPS, MATSOLVERSPOOLES, PCFactorSetMatSolverPackage(), MatSolverPackage
M*/
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "MatGetFactor_seqaij_umfpack"
PetscErrorCode MatGetFactor_seqaij_umfpack(Mat A,MatFactorType ftype,Mat *F)
{
Mat B;
Mat_UMFPACK *lu;
PetscErrorCode ierr;
PetscInt m=A->rmap->n,n=A->cmap->n,idx;
const char *strategy[]={"AUTO","UNSYMMETRIC","SYMMETRIC"},
*scale[]={"NONE","SUM","MAX"};
PetscBool flg;
PetscFunctionBegin;
/* Create the factorization matrix F */
ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr);
ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr);
ierr = MatSetType(B,((PetscObject)A)->type_name);CHKERRQ(ierr);
ierr = MatSeqAIJSetPreallocation(B,0,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscNewLog(B,Mat_UMFPACK,&lu);CHKERRQ(ierr);
B->spptr = lu;
B->ops->lufactorsymbolic = MatLUFactorSymbolic_UMFPACK;
B->ops->destroy = MatDestroy_UMFPACK;
B->ops->view = MatView_UMFPACK;
ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatFactorGetSolverPackage_C","MatFactorGetSolverPackage_seqaij_umfpack",MatFactorGetSolverPackage_seqaij_umfpack);CHKERRQ(ierr);
B->factortype = MAT_FACTOR_LU;
B->assembled = PETSC_TRUE; /* required by -ksp_view */
B->preallocated = PETSC_TRUE;
/* initializations */
/* ------------------------------------------------*/
/* get the default control parameters */
umfpack_UMF_defaults(lu->Control);
lu->perm_c = PETSC_NULL; /* use defaul UMFPACK col permutation */
lu->Control[UMFPACK_IRSTEP] = 0; /* max num of iterative refinement steps to attempt */
ierr = PetscOptionsBegin(((PetscObject)A)->comm,((PetscObject)A)->prefix,"UMFPACK Options","Mat");CHKERRQ(ierr);
/* Control parameters used by reporting routiones */
ierr = PetscOptionsReal("-mat_umfpack_prl","Control[UMFPACK_PRL]","None",lu->Control[UMFPACK_PRL],&lu->Control[UMFPACK_PRL],PETSC_NULL);CHKERRQ(ierr);
/* Control parameters for symbolic factorization */
ierr = PetscOptionsEList("-mat_umfpack_strategy","ordering and pivoting strategy","None",strategy,3,strategy[0],&idx,&flg);CHKERRQ(ierr);
if (flg) {
switch (idx){
case 0: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_AUTO; break;
case 1: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_UNSYMMETRIC; break;
case 2: lu->Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC; break;
}
}
ierr = PetscOptionsEList("-mat_umfpack_ordering","Internal ordering method","None",UmfpackOrderingTypes,sizeof UmfpackOrderingTypes/sizeof UmfpackOrderingTypes[0],UmfpackOrderingTypes[(int)lu->Control[UMFPACK_ORDERING]],&idx,&flg);CHKERRQ(ierr);
if (flg) lu->Control[UMFPACK_ORDERING] = (int)idx;
ierr = PetscOptionsReal("-mat_umfpack_dense_col","Control[UMFPACK_DENSE_COL]","None",lu->Control[UMFPACK_DENSE_COL],&lu->Control[UMFPACK_DENSE_COL],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_dense_row","Control[UMFPACK_DENSE_ROW]","None",lu->Control[UMFPACK_DENSE_ROW],&lu->Control[UMFPACK_DENSE_ROW],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_amd_dense","Control[UMFPACK_AMD_DENSE]","None",lu->Control[UMFPACK_AMD_DENSE],&lu->Control[UMFPACK_AMD_DENSE],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_block_size","Control[UMFPACK_BLOCK_SIZE]","None",lu->Control[UMFPACK_BLOCK_SIZE],&lu->Control[UMFPACK_BLOCK_SIZE],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_fixq","Control[UMFPACK_FIXQ]","None",lu->Control[UMFPACK_FIXQ],&lu->Control[UMFPACK_FIXQ],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_aggressive","Control[UMFPACK_AGGRESSIVE]","None",lu->Control[UMFPACK_AGGRESSIVE],&lu->Control[UMFPACK_AGGRESSIVE],PETSC_NULL);CHKERRQ(ierr);
/* Control parameters used by numeric factorization */
ierr = PetscOptionsReal("-mat_umfpack_pivot_tolerance","Control[UMFPACK_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_PIVOT_TOLERANCE],&lu->Control[UMFPACK_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-mat_umfpack_sym_pivot_tolerance","Control[UMFPACK_SYM_PIVOT_TOLERANCE]","None",lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],&lu->Control[UMFPACK_SYM_PIVOT_TOLERANCE],PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsEList("-mat_umfpack_scale","Control[UMFPACK_SCALE]","None",scale,3,scale[0],&idx,&flg);CHKERRQ(ierr);
if (flg) {
switch (idx){
case 0: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_NONE; break;
case 1: lu->Control[UMFPACK_SCALE] = UMFPACK_SCALE_SUM; break;
//.........这里部分代码省略.........
开发者ID:Kun-Qu，项目名称:petsc，代码行数:101，代码来源:umfpack.c
-
六六分期app的软件客服如何联系？不知道吗？加qq群【895510560】即可！标题：六六分期
阅读：18213|2023-10-27
-
今天小编告诉大家如何处理win10系统火狐flash插件总是崩溃的问题，可能很多用户都不知
阅读：9656|2022-11-06
-
今天小编告诉大家如何对win10系统删除桌面回收站图标进行设置，可能很多用户都不知道
阅读：8168|2022-11-06
-
今天小编告诉大家如何对win10系统电脑设置节能降温的设置方法，想必大家都遇到过需要
阅读：8543|2022-11-06
-
我们在使用xp系统的过程中,经常需要对xp系统无线网络安装向导设置进行设置，可能很多
阅读：8449|2022-11-06
-
今天小编告诉大家如何处理win7系统玩cf老是与主机连接不稳定的问题，可能很多用户都不
阅读：9375|2022-11-06
-
电脑对日常生活的重要性小编就不多说了，可是一旦碰到win7系统设置cf烟雾头的问题，很
阅读：8418|2022-11-06
-
我们在日常使用电脑的时候，有的小伙伴们可能在打开应用的时候会遇见提示应用程序无法
阅读：7855|2022-11-06
-
今天小编告诉大家如何对win7系统打开vcf文件进行设置，可能很多用户都不知道怎么对win
阅读：8403|2022-11-06
-
今天小编告诉大家如何对win10系统s4开启USB调试模式进行设置，可能很多用户都不知道怎
阅读：7391|2022-11-06
|
客服电话
APP下载
官方微信
请发表评论