当前位置: 首页>>代码示例>>C++>>正文


C++ KSPCreate函数代码示例

本文整理汇总了C++中KSPCreate函数的典型用法代码示例。如果您正苦于以下问题:C++ KSPCreate函数的具体用法?C++ KSPCreate怎么用?C++ KSPCreate使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。


在下文中一共展示了KSPCreate函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。

示例1: petsc_solve

int petsc_solve(int n_, double complex *A_, double complex *b_, double complex *x_)
{
    Vec            x, b;
    Mat            A;
    KSP            ksp;
    PC             pc;
    PetscReal      norm, tol=1.e-14;
    PetscErrorCode ierr;
    PetscInt       i, j, n = n_, col[n_], its;
    PetscScalar    neg_one = -1.0, one = 1.0, value[n_], *x_array;

    ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject) x, "Solution");CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
    ierr = VecSetFromOptions(x);CHKERRQ(ierr);
    ierr = VecDuplicate(x,&b);CHKERRQ(ierr);

    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
    ierr = MatSetFromOptions(A);CHKERRQ(ierr);
    ierr = MatSetOption(A,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatSetUp(A);CHKERRQ(ierr);

    for (i=0; i<n; i++) col[i] = i;
    printf("  Converting matrix to PETSc\n");
    ierr = MatSetValues(A,n,col,n,col,A_,INSERT_VALUES);CHKERRQ(ierr);
    printf("    Done\n");
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

    for (i=0; i<n; i++) value[i] = b_[i];
    ierr = VecSetValues(b, n, col, value, INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

    ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

    ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
    ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
    // For a full list, see:
    // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPType.html
    ierr = KSPSetType(ksp, KSPCGS);CHKERRQ(ierr);

    printf("  Solving...\n");
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
    printf("  Done\n");

//    ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

    ierr = VecGetArray(x, &x_array);CHKERRQ(ierr);
    for (i=0; i<n; i++) x_[i] = x_array[i];

    ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Iterations %D\n", its);CHKERRQ(ierr);

    ierr = VecDestroy(&x);CHKERRQ(ierr);
    ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
    ierr = KSPDestroy(&ksp);CHKERRQ(ierr);

    return 0;
}
开发者ID:certik,项目名称:hfsolver,代码行数:66,代码来源:c_petsc.c

示例2: TaylorGalerkinStepIIMassEnergy

/* Notice that this requires the previous momentum solution.

The element stiffness matrix for the identity in linear elements is

  1  /2 1 1\
  -  |1 2 1|
  12 \1 1 2/

  no matter what the shape of the triangle. */
PetscErrorCode TaylorGalerkinStepIIMassEnergy(DM da, UserContext *user)
{
  MPI_Comm       comm;
  Mat            mat;
  Vec            rhs_m, rhs_e;
  PetscScalar    identity[9] = {0.16666666667, 0.08333333333, 0.08333333333,
                                0.08333333333, 0.16666666667, 0.08333333333,
                                0.08333333333, 0.08333333333, 0.16666666667};
  PetscScalar    *u_n,       *v_n,     *p_n,     *t_n,     *mu_n,    *kappa_n;
  PetscScalar    *rho_n,     *rho_u_n, *rho_v_n, *rho_e_n;
  PetscScalar    *u_phi,     *v_phi;
  PetscScalar    *rho_u_np1, *rho_v_np1;
  PetscInt       idx[3];
  PetscScalar    psi_x[3], psi_y[3];
  PetscScalar    values_m[3];
  PetscScalar    values_e[3];
  PetscScalar    phi = user->phi;
  PetscScalar    mu, kappa, tau_xx, tau_xy, tau_yy, q_x, q_y;
  PetscReal      hx, hy, area;
  KSP            ksp;
  const PetscInt *necon;
  PetscInt       j, k, e, ne, nc, mx, my;
  PetscErrorCode ierr;

  PetscFunctionBeginUser;
  ierr = PetscObjectGetComm((PetscObject) da, &comm);CHKERRQ(ierr);
  ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
  ierr = DMCreateMatrix(da, &mat);CHKERRQ(ierr);
  ierr = MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &rhs_m);CHKERRQ(ierr);
  ierr = DMGetGlobalVector(da, &rhs_e);CHKERRQ(ierr);
  ierr = KSPCreate(comm, &ksp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  ierr = DMDAGetInfo(da, 0, &mx, &my, 0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
  hx   = 1.0 / (PetscReal)(mx-1);
  hy   = 1.0 / (PetscReal)(my-1);
  area = 0.5*hx*hy;
  ierr = VecGetArray(user->sol_n.u,       &u_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.v,       &v_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.p,       &p_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.t,       &t_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->mu,            &mu_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->kappa,         &kappa_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.rho,     &rho_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.rho_u,   &rho_u_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.rho_v,   &rho_v_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_n.rho_e,   &rho_e_n);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.u,     &u_phi);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_phi.v,     &v_phi);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_np1.rho_u, &rho_u_np1);CHKERRQ(ierr);
  ierr = VecGetArray(user->sol_np1.rho_v, &rho_v_np1);CHKERRQ(ierr);
  ierr = DMDAGetElements(da, &ne, &nc, &necon);CHKERRQ(ierr);
  for (e = 0; e < ne; e++) {
    for (j = 0; j < 3; j++) {
      idx[j]      = necon[3*e+j];
      values_m[j] = 0.0;
      values_e[j] = 0.0;
    }
    /* Get basis function deriatives (we need the orientation of the element here) */
    if (idx[1] > idx[0]) {
      psi_x[0] = -hy; psi_x[1] =  hy; psi_x[2] = 0.0;
      psi_y[0] = -hx; psi_y[1] = 0.0; psi_y[2] =  hx;
    } else {
      psi_x[0] =  hy; psi_x[1] = -hy; psi_x[2] = 0.0;
      psi_y[0] =  hx; psi_y[1] = 0.0; psi_y[2] = -hx;
    }
    /*  <\nabla\psi, F^*>: Divergence of the predicted convective fluxes */
    for (j = 0; j < 3; j++) {
      values_m[j] += (psi_x[j]*(phi*rho_u_np1[idx[j]] + rho_u_n[idx[j]]) + psi_y[j]*(rho_v_np1[idx[j]] + rho_v_n[idx[j]]))/3.0;
      values_e[j] += values_m[j]*((rho_e_n[idx[j]] + p_n[idx[j]]) / rho_n[idx[j]]);
    }
    /*  -<\nabla\psi, F^n_v>: Divergence of the viscous fluxes */
    for (j = 0; j < 3; j++) {
      /* \tau_{xx} = 2/3 \mu(T) (2 {\partial u\over\partial x} - {\partial v\over\partial y}) */
      /* \tau_{xy} =     \mu(T) (  {\partial u\over\partial y} + {\partial v\over\partial x}) */
      /* \tau_{yy} = 2/3 \mu(T) (2 {\partial v\over\partial y} - {\partial u\over\partial x}) */
      /* q_x       = -\kappa(T) {\partial T\over\partial x} */
      /* q_y       = -\kappa(T) {\partial T\over\partial y} */

      /* above code commeted out - causing ininitialized variables. */
      q_x =0; q_y =0;

      mu     = 0.0;
      kappa  = 0.0;
      tau_xx = 0.0;
      tau_xy = 0.0;
      tau_yy = 0.0;
      for (k = 0; k < 3; k++) {
        mu     += mu_n[idx[k]];
        kappa  += kappa_n[idx[k]];
//.........这里部分代码省略.........
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:101,代码来源:ex31.c

示例3: PCISSetUp


//.........这里部分代码省略.........
  ierr = VecScatterCreate(pcis->vec1_N,pcis->is_B_local,pcis->vec1_B,(IS)0,&pcis->N_to_B);CHKERRQ(ierr);
  ierr = VecScatterCreate(pcis->vec1_global,pcis->is_B_global,pcis->vec1_B,(IS)0,&pcis->global_to_B);CHKERRQ(ierr);

  /* Creating scaling "matrix" D */
  ierr = PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_is_use_stiffness_scaling",&pcis->use_stiffness_scaling,NULL);CHKERRQ(ierr);
  if (!pcis->D) {
    ierr = VecDuplicate(pcis->vec1_B,&pcis->D);CHKERRQ(ierr);
    if (!pcis->use_stiffness_scaling) {
      ierr = VecSet(pcis->D,pcis->scaling_factor);CHKERRQ(ierr);
    } else {
      ierr = MatGetDiagonal(matis->A,pcis->vec1_N);CHKERRQ(ierr);
      ierr = VecScatterBegin(pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
      ierr = VecScatterEnd  (pcis->N_to_B,pcis->vec1_N,pcis->D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    }
  }
  ierr = VecCopy(pcis->D,pcis->vec1_B);CHKERRQ(ierr);
  ierr = MatCreateVecs(pc->pmat,&counter,0);CHKERRQ(ierr); /* temporary auxiliar vector */
  ierr = VecSet(counter,0.0);CHKERRQ(ierr);
  ierr = VecScatterBegin(pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->global_to_B,pcis->vec1_B,counter,ADD_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
  ierr = VecScatterBegin(pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecScatterEnd  (pcis->global_to_B,counter,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
  ierr = VecPointwiseDivide(pcis->D,pcis->D,pcis->vec1_B);CHKERRQ(ierr);
  ierr = VecDestroy(&counter);CHKERRQ(ierr);

  /* See historical note 01, at the bottom of this file. */

  /*
    Creating the KSP contexts for the local Dirichlet and Neumann problems.
  */
  if (pcis->computesolvers) {
    PC pc_ctx;
    /* Dirichlet */
    ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_D);CHKERRQ(ierr);
    ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr);
    ierr = KSPSetOperators(pcis->ksp_D,pcis->A_II,pcis->A_II);CHKERRQ(ierr);
    ierr = KSPSetOptionsPrefix(pcis->ksp_D,"is_localD_");CHKERRQ(ierr);
    ierr = KSPGetPC(pcis->ksp_D,&pc_ctx);CHKERRQ(ierr);
    ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
    ierr = KSPSetType(pcis->ksp_D,KSPPREONLY);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(pcis->ksp_D);CHKERRQ(ierr);
    /* the vectors in the following line are dummy arguments, just telling the KSP the vector size. Values are not used */
    ierr = KSPSetUp(pcis->ksp_D);CHKERRQ(ierr);
    /* Neumann */
    ierr = KSPCreate(PETSC_COMM_SELF,&pcis->ksp_N);CHKERRQ(ierr);
    ierr = PetscObjectIncrementTabLevel((PetscObject)pcis->ksp_N,(PetscObject)pc,1);CHKERRQ(ierr);
    ierr = KSPSetOperators(pcis->ksp_N,matis->A,matis->A);CHKERRQ(ierr);
    ierr = KSPSetOptionsPrefix(pcis->ksp_N,"is_localN_");CHKERRQ(ierr);
    ierr = KSPGetPC(pcis->ksp_N,&pc_ctx);CHKERRQ(ierr);
    ierr = PCSetType(pc_ctx,PCLU);CHKERRQ(ierr);
    ierr = KSPSetType(pcis->ksp_N,KSPPREONLY);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(pcis->ksp_N);CHKERRQ(ierr);
    {
      PetscBool damp_fixed                    = PETSC_FALSE,
                remove_nullspace_fixed        = PETSC_FALSE,
                set_damping_factor_floating   = PETSC_FALSE,
                not_damp_floating             = PETSC_FALSE,
                not_remove_nullspace_floating = PETSC_FALSE;
      PetscReal fixed_factor,
                floating_factor;

      ierr = PetscOptionsGetReal(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&fixed_factor,&damp_fixed);CHKERRQ(ierr);
      if (!damp_fixed) fixed_factor = 0.0;
      ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_damp_fixed",&damp_fixed,NULL);CHKERRQ(ierr);

      ierr = PetscOptionsGetBool(((PetscObject)pc_ctx)->prefix,"-pc_is_remove_nullspace_fixed",&remove_nullspace_fixed,NULL);CHKERRQ(ierr);
开发者ID:fengyuqi,项目名称:petsc,代码行数:67,代码来源:pcis.c

示例4: main

int main(int argc, char **argv)
{
  PetscErrorCode ierr;
  PetscInt       n=10000,its,dfid=1;
  Vec            x,b,u;
  Mat            A;
  KSP            ksp;
  PC             pc,pcnoise;
  PCNoise_Ctx    ctx={0,NULL};
  PetscReal      eta=0.1,norm;
  PetscScalar(*diagfunc)(PetscInt,PetscInt);

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);CHKERRQ(ierr);

  /* Process command line options */
  ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetReal(NULL,"-eta",&eta,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,"-diagfunc",&dfid,NULL);CHKERRQ(ierr);
  switch(dfid){
    case 1:
      diagfunc = diagFunc1;
      break;
    case 2:
      diagfunc = diagFunc2;
      break;
    case 3:
      diagfunc = diagFunc3;
      break;
    default:
      SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Unrecognized diagfunc option");
  }

  /* Create a diagonal matrix with a given distribution of diagonal elements */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatSetUp(A);CHKERRQ(ierr);
  ierr = AssembleDiagonalMatrix(A,diagfunc);CHKERRQ(ierr);

  /* Allocate vectors and manufacture an exact solution and rhs */
  ierr = MatCreateVecs(A,&x,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)x,"Computed Solution");CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&b,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)b,"RHS");CHKERRQ(ierr);
  ierr = MatCreateVecs(A,&u,NULL);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)u,"Reference Solution");CHKERRQ(ierr);
  ierr = VecSet(u,1.0);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /* Create a KSP object */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

  /* Set up a composite preconditioner */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCCOMPOSITE);CHKERRQ(ierr); /* default composite with single Identity PC */
  ierr = PCCompositeSetType(pc,PC_COMPOSITE_ADDITIVE);CHKERRQ(ierr);
  ierr = PCCompositeAddPC(pc,PCNONE);CHKERRQ(ierr);
  if(eta > 0){
    ierr = PCCompositeAddPC(pc,PCSHELL);CHKERRQ(ierr);
    ierr = PCCompositeGetPC(pc,1,&pcnoise);CHKERRQ(ierr);
    ctx.eta = eta;
    ierr = PCShellSetContext(pcnoise,&ctx);CHKERRQ(ierr);
    ierr = PCShellSetApply(pcnoise,PCApply_Noise);CHKERRQ(ierr);
    ierr = PCShellSetSetUp(pcnoise,PCSetup_Noise);CHKERRQ(ierr);
    ierr = PCShellSetDestroy(pcnoise,PCDestroy_Noise);CHKERRQ(ierr);
    ierr = PCShellSetName(pcnoise,"Noise PC");CHKERRQ(ierr);
  }

  /* Set KSP from options (this can override the PC just defined) */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

  /* Solve */
  ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);

  /* Compute error */
  ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
  ierr = PetscObjectSetName((PetscObject)x,"Error");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",(double)norm,its);CHKERRQ(ierr);

  /* Destroy objects and finalize */
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&x);CHKERRQ(ierr);
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&u);CHKERRQ(ierr);
  PetscFinalize();

  return 0;
}
开发者ID:pombredanne,项目名称:petsc,代码行数:92,代码来源:ex60.c

示例5: 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 */
  KSP            *subksp;     /* array of local KSP contexts on this processor */
  PC             pc;           /* PC context */
  PC             subpc;        /* PC context for subdomain */
  PetscReal      norm;         /* norm of solution error */
  PetscErrorCode ierr;
  PetscInt       i,j,Ii,J,*blks,m = 8,n;
  PetscMPIInt    rank,size;
  PetscInt       its,nlocal,first,Istart,Iend;
  PetscScalar    v,one = 1.0,none = -1.0;
  PetscBool      isbjacobi,flg = PETSC_FALSE;

  PetscInitialize(&argc,&args,(char*)0,help);
  ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
  n    = m+2;

  /* -------------------------------------------------------------------
         Compute the matrix and right-hand-side vector that define
         the linear system, Ax = b.
     ------------------------------------------------------------------- */

  /*
     Create and assemble parallel matrix
  */
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr);
  ierr = MatSeqAIJSetPreallocation(A,5,NULL);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; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);}
    if (i<m-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; 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);

  /*
     Create parallel vectors
  */
  ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
  ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
  ierr = VecSetFromOptions(u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);

  /*
     Set exact solution; then compute right-hand-side vector.
  */
  ierr = VecSet(u,one);CHKERRQ(ierr);
  ierr = MatMult(A,u,b);CHKERRQ(ierr);

  /*
     Create linear solver context
  */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /*
     Set operators. Here the matrix that defines the linear system
     also serves as the preconditioning matrix.
  */
  ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

  /*
     Set default preconditioner for this program to be block Jacobi.
     This choice can be overridden at runtime with the option
        -pc_type <type>
  */
  ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
  ierr = PCSetType(pc,PCBJACOBI);CHKERRQ(ierr);


  /* -------------------------------------------------------------------
                   Define the problem decomposition
     ------------------------------------------------------------------- */

  /*
     Call PCBJacobiSetTotalBlocks() to set individually the size of
     each block in the preconditioner.  This could also be done with
     the runtime option
         -pc_bjacobi_blocks <blocks>
     Also, see the command PCBJacobiSetLocalBlocks() to set the
     local blocks.

      Note: The default decomposition is 1 block per processor.
  */
  ierr = PetscMalloc(m*sizeof(PetscInt),&blks);CHKERRQ(ierr);
  for (i=0; i<m; i++) blks[i] = n;
  ierr = PCBJacobiSetTotalBlocks(pc,m,blks);CHKERRQ(ierr);
  ierr = PetscFree(blks);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:feelpp,项目名称:debian-petsc,代码行数:101,代码来源:ex7.c

示例6: PCBDDCNullSpaceAssembleCorrection


//.........这里部分代码省略.........
      ierr = VecAssemblyBegin(shell_ctx->work_small_1);CHKERRQ(ierr);
      ierr = VecAssemblyEnd(shell_ctx->work_small_1);CHKERRQ(ierr);
      ierr = MatSolve(small_mat,shell_ctx->work_small_1,shell_ctx->work_small_2);CHKERRQ(ierr);
      ierr = VecGetArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
      for (i=0;i<basis_size;i++) {
        array_mat[i*basis_size+k]=array[i];
      }
      ierr = VecRestoreArrayRead(shell_ctx->work_small_2,(const PetscScalar**)&array);CHKERRQ(ierr);
    }
    ierr = MatCreateSeqDense(PETSC_COMM_SELF,basis_size,basis_size,array_mat,&inv_small_mat);CHKERRQ(ierr);
    ierr = MatMatMult(shell_ctx->basis_mat,inv_small_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&shell_ctx->Lbasis_mat);CHKERRQ(ierr);
    ierr = PetscFree(array_mat);CHKERRQ(ierr);
    ierr = MatDestroy(&inv_small_mat);CHKERRQ(ierr);
    ierr = MatDestroy(&small_mat);CHKERRQ(ierr);
    ierr = MatScale(shell_ctx->Kbasis_mat,m_one);CHKERRQ(ierr);

    /* Rebuild local PC */
    ierr = KSPGetPC(local_ksp,&shell_ctx->local_pc);CHKERRQ(ierr);
    ierr = PetscObjectReference((PetscObject)shell_ctx->local_pc);CHKERRQ(ierr);
    ierr = PCCreate(PETSC_COMM_SELF,&newpc);CHKERRQ(ierr);
    ierr = PCSetOperators(newpc,local_mat,local_mat);CHKERRQ(ierr);
    ierr = PCSetType(newpc,PCSHELL);CHKERRQ(ierr);
    ierr = PCShellSetContext(newpc,shell_ctx);CHKERRQ(ierr);
    ierr = PCShellSetApply(newpc,PCBDDCApplyNullSpaceCorrectionPC);CHKERRQ(ierr);
    ierr = PCShellSetDestroy(newpc,PCBDDCDestroyNullSpaceCorrectionPC);CHKERRQ(ierr);
    ierr = PCSetUp(newpc);CHKERRQ(ierr);
    ierr = KSPSetPC(local_ksp,newpc);CHKERRQ(ierr);
    ierr = PCDestroy(&newpc);CHKERRQ(ierr);
    ierr = KSPSetUp(local_ksp);CHKERRQ(ierr);
  }
  /* test */
  if (pcbddc->dbg_flag && basis_dofs) {
    KSP         check_ksp;
    PC          check_pc;
    Mat         test_mat;
    Vec         work3;
    PetscReal   test_err,lambda_min,lambda_max;
    PetscBool   setsym,issym=PETSC_FALSE;
    PetscInt    tabs;

    ierr = PetscViewerASCIIGetTab(pcbddc->dbg_viewer,&tabs);CHKERRQ(ierr);
    ierr = KSPGetPC(local_ksp,&check_pc);CHKERRQ(ierr);
    ierr = VecDuplicate(shell_ctx->work_full_1,&work1);CHKERRQ(ierr);
    ierr = VecDuplicate(shell_ctx->work_full_1,&work2);CHKERRQ(ierr);
    ierr = VecDuplicate(shell_ctx->work_full_1,&work3);CHKERRQ(ierr);
    ierr = VecSetRandom(shell_ctx->work_small_1,NULL);CHKERRQ(ierr);
    ierr = MatMult(shell_ctx->basis_mat,shell_ctx->work_small_1,work1);CHKERRQ(ierr);
    ierr = VecCopy(work1,work2);CHKERRQ(ierr);
    ierr = MatMult(local_mat,work1,work3);CHKERRQ(ierr);
    ierr = PCApply(check_pc,work3,work1);CHKERRQ(ierr);
    ierr = VecAXPY(work1,m_one,work2);CHKERRQ(ierr);
    ierr = VecNorm(work1,NORM_INFINITY,&test_err);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace correction for ",PetscGlobalRank);CHKERRQ(ierr);
    ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_FALSE);CHKERRQ(ierr);
    if (isdir) {
      ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Dirichlet ");CHKERRQ(ierr);
    } else {
      ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Neumann ");CHKERRQ(ierr);
    }
    ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"solver is :%1.14e\n",test_err);CHKERRQ(ierr);
    ierr = PetscViewerASCIISetTab(pcbddc->dbg_viewer,tabs);CHKERRQ(ierr);
    ierr = PetscViewerASCIIUseTabs(pcbddc->dbg_viewer,PETSC_TRUE);CHKERRQ(ierr);

    ierr = MatTransposeMatMult(shell_ctx->Lbasis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&test_mat);CHKERRQ(ierr);
    ierr = MatShift(test_mat,one);CHKERRQ(ierr);
    ierr = MatNorm(test_mat,NORM_INFINITY,&test_err);CHKERRQ(ierr);
    ierr = MatDestroy(&test_mat);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for nullspace matrices is :%1.14e\n",PetscGlobalRank,test_err);CHKERRQ(ierr);

    /* Create ksp object suitable for extreme eigenvalues' estimation */
    ierr = KSPCreate(PETSC_COMM_SELF,&check_ksp);CHKERRQ(ierr);
    ierr = KSPSetErrorIfNotConverged(check_ksp,pc->erroriffailure);CHKERRQ(ierr);
    ierr = KSPSetOperators(check_ksp,local_mat,local_mat);CHKERRQ(ierr);
    ierr = KSPSetTolerances(check_ksp,1.e-8,1.e-8,PETSC_DEFAULT,basis_dofs);CHKERRQ(ierr);
    ierr = KSPSetComputeSingularValues(check_ksp,PETSC_TRUE);CHKERRQ(ierr);
    ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr);
    if (issym) {
      ierr = KSPSetType(check_ksp,KSPCG);CHKERRQ(ierr);
    }
    ierr = KSPSetPC(check_ksp,check_pc);CHKERRQ(ierr);
    ierr = KSPSetUp(check_ksp);CHKERRQ(ierr);
    ierr = VecSetRandom(work1,NULL);CHKERRQ(ierr);
    ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
    ierr = KSPSolve(check_ksp,work2,work2);CHKERRQ(ierr);
    ierr = VecAXPY(work2,m_one,work1);CHKERRQ(ierr);
    ierr = VecNorm(work2,NORM_INFINITY,&test_err);CHKERRQ(ierr);
    ierr = KSPComputeExtremeSingularValues(check_ksp,&lambda_max,&lambda_min);CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(check_ksp,&k);CHKERRQ(ierr);
    ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d error for adapted KSP %1.14e (it %d, eigs %1.6e %1.6e)\n",PetscGlobalRank,test_err,k,lambda_min,lambda_max);CHKERRQ(ierr);
    ierr = KSPDestroy(&check_ksp);CHKERRQ(ierr);
    ierr = VecDestroy(&work1);CHKERRQ(ierr);
    ierr = VecDestroy(&work2);CHKERRQ(ierr);
    ierr = VecDestroy(&work3);CHKERRQ(ierr);
  }
  /* all processes shoud call this, even the void ones */
  if (pcbddc->dbg_flag) {
    ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}
开发者ID:plguhur,项目名称:petsc,代码行数:101,代码来源:bddcnullspace.c

示例7: main


//.........这里部分代码省略.........
      ierr    = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
      sxy.c   = 1; sxy_m.c = 0; /* spin 1, 0 */
      val     = -PETSC_i*uxy2; valconj = PetscConj(val);
      ierr    = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
      ierr    = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
      sxy.c   = 1; sxy_m.c = 1; /* spin 1, 1 */
      val     = PetscConj(uxy2); valconj = PetscConj(val);
      ierr    = MatSetValuesStencil(H, 1, &sxy_m, 1, &sxy, &val, ADD_VALUES);CHKERRQ(ierr);
      ierr    = MatSetValuesStencil(H, 1, &sxy, 1, &sxy_m, &valconj, ADD_VALUES);CHKERRQ(ierr);
    }
  }

  ierr = MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  ierr = MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

  /* scale H */
  ierr = MatScale(H, 1./(2.*h));CHKERRQ(ierr);

  /* it looks like H is Hermetian */
  /* construct normal equations */
  ierr = MatMatMult(H, H, MAT_INITIAL_MATRIX, 1., &HtH);CHKERRQ(ierr);

  /* permutation matrix to check whether H and HtH are identical to the ones in the paper */
/*   Mat perm; */
/*   ierr = DMCreateMatrix(da, &perm);CHKERRQ(ierr); */
/*   PetscInt row, col; */
/*   PetscScalar one = 1.0; */
/*   for (PetscInt i=0; i<n; i++) { */
/*     for (PetscInt j=0; j<n; j++) { */
/*       row = (i*n+j)*2; col = i*n+j; */
/*       ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
/*       row = (i*n+j)*2+1; col = i*n+j + n*n; */
/*       ierr = MatSetValues(perm, 1, &row, 1, &col, &one, INSERT_VALUES);CHKERRQ(ierr); */
/*     } */
/*   } */
/*   ierr = MatAssemblyBegin(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */
/*   ierr = MatAssemblyEnd(perm, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); */

/*   Mat Hperm; */
/*   ierr = MatPtAP(H, perm, MAT_INITIAL_MATRIX, 1.0, &Hperm);CHKERRQ(ierr); */
/*   ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix H after construction\n");CHKERRQ(ierr); */
/*   ierr = MatView(Hperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */

/*   Mat HtHperm; */
/*   ierr = MatPtAP(HtH, perm, MAT_INITIAL_MATRIX, 1.0, &HtHperm);CHKERRQ(ierr); */
/*   ierr = PetscPrintf(PETSC_COMM_WORLD, "Matrix HtH:\n");CHKERRQ(ierr); */
/*   ierr = MatView(HtHperm, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */

  /* right hand side */
  ierr = DMCreateGlobalVector(da, &b);CHKERRQ(ierr);
  ierr = VecSet(b,0.0);CHKERRQ(ierr);
  ierr = VecSetValues(b, 1, ix, vals, INSERT_VALUES);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
/*   ierr = VecSetRandom(b, rctx);CHKERRQ(ierr); */
  ierr = VecDuplicate(b, &Htb);CHKERRQ(ierr);
  ierr = MatMultTranspose(H, b, Htb);CHKERRQ(ierr);

  /* construct solver */
  ierr = KSPCreate(PETSC_COMM_WORLD,&kspmg);CHKERRQ(ierr);
  ierr = KSPSetType(kspmg, KSPCG);CHKERRQ(ierr);

  ierr = KSPGetPC(kspmg,&pcmg);CHKERRQ(ierr);
  ierr = PCSetType(pcmg,PCASA);CHKERRQ(ierr);

  /* maybe user wants to override some of the choices */
  ierr = KSPSetFromOptions(kspmg);CHKERRQ(ierr);

  ierr = KSPSetOperators(kspmg, HtH, HtH, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);

  ierr = DMDASetRefinementFactor(da, 3, 3, 3);CHKERRQ(ierr);
  ierr = PCSetDM(pcmg,da);CHKERRQ(ierr);

  ierr = PCASASetTolerances(pcmg, 1.e-6, 1.e-10,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr);

  ierr = VecDuplicate(b, &xvec);CHKERRQ(ierr);
  ierr = VecSet(xvec, 0.0);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                      Solve the linear system
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  ierr = KSPSolve(kspmg, Htb, xvec);CHKERRQ(ierr);

/*   ierr = VecView(xvec, PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD));CHKERRQ(ierr); */

  ierr = KSPDestroy(&kspmg);CHKERRQ(ierr);
  ierr = VecDestroy(&xvec);CHKERRQ(ierr);

  /*   seems to be destroyed by KSPDestroy */
  ierr = VecDestroy(&b);CHKERRQ(ierr);
  ierr = VecDestroy(&Htb);CHKERRQ(ierr);
  ierr = MatDestroy(&HtH);CHKERRQ(ierr);
  ierr = MatDestroy(&H);CHKERRQ(ierr);

  ierr = DMDestroy(&da);CHKERRQ(ierr);
  ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
  ierr = PetscFinalize();
  return 0;
}
开发者ID:hsahasra,项目名称:petsc-magma-dense-mat,代码行数:101,代码来源:ex39.c

示例8: main

int main(int argc,char **argv)
{
  PetscErrorCode ierr;
  PetscInt       its,n,Nx=PETSC_DECIDE,Ny=PETSC_DECIDE,nlocal;
  PetscMPIInt    size;
  PetscScalar    one = 1.0;
  PetscInt       mx,my;
  Mat            A;
  GridCtx        fine_ctx;
  KSP            ksp;
  PetscBool      flg;

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  /* set up discretization matrix for fine grid */
  fine_ctx.mx = 9; fine_ctx.my = 9;
  ierr        = PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);CHKERRQ(ierr);
  if (flg) fine_ctx.mx = mx;
  ierr = PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);CHKERRQ(ierr);
  if (flg) fine_ctx.my = my;
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Fine grid size %D by %D\n",fine_ctx.mx,fine_ctx.my);CHKERRQ(ierr);
  n    = fine_ctx.mx*fine_ctx.my;

  MPI_Comm_size(PETSC_COMM_WORLD,&size);
  ierr = PetscOptionsGetInt(NULL,NULL,"-Nx",&Nx,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,NULL,"-Ny",&Ny,NULL);CHKERRQ(ierr);

  /* Set up distributed array for fine grid */
  ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,fine_ctx.mx,
                      fine_ctx.my,Nx,Ny,1,1,NULL,NULL,&fine_ctx.da);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(fine_ctx.da,&fine_ctx.x);CHKERRQ(ierr);
  ierr = VecDuplicate(fine_ctx.x,&fine_ctx.b);CHKERRQ(ierr);
  ierr = VecGetLocalSize(fine_ctx.x,&nlocal);CHKERRQ(ierr);
  ierr = DMCreateLocalVector(fine_ctx.da,&fine_ctx.localX);CHKERRQ(ierr);
  ierr = VecDuplicate(fine_ctx.localX,&fine_ctx.localF);CHKERRQ(ierr);
  ierr = MatCreateAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,NULL,3,NULL,&A);CHKERRQ(ierr);
  ierr = FormJacobian_Grid(&fine_ctx,&A);CHKERRQ(ierr);

  /* create linear solver */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

  /* set values for rhs vector */
  ierr = VecSet(fine_ctx.b,one);CHKERRQ(ierr);

  /* set options, then solve system */
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); /* calls PCSetFromOptions_ML if 'pc_type=ml' */
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,fine_ctx.b,fine_ctx.x);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %D\n",its);CHKERRQ(ierr);

  /* free data structures */
  ierr = VecDestroy(&fine_ctx.x);CHKERRQ(ierr);
  ierr = VecDestroy(&fine_ctx.b);CHKERRQ(ierr);
  ierr = DMDestroy(&fine_ctx.da);CHKERRQ(ierr);
  ierr = VecDestroy(&fine_ctx.localX);CHKERRQ(ierr);
  ierr = VecDestroy(&fine_ctx.localF);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);

  ierr = PetscFinalize();
  return ierr;
}
开发者ID:tom-klotz,项目名称:petsc,代码行数:62,代码来源:ex26.c

示例9: TaoCreate_TRON

PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
{
  TAO_TRON       *tron;
  PetscErrorCode ierr;
  const char     *morethuente_type = TAOLINESEARCHMT;

  PetscFunctionBegin;
  tao->ops->setup = TaoSetup_TRON;
  tao->ops->solve = TaoSolve_TRON;
  tao->ops->view = TaoView_TRON;
  tao->ops->setfromoptions = TaoSetFromOptions_TRON;
  tao->ops->destroy = TaoDestroy_TRON;
  tao->ops->computedual = TaoComputeDual_TRON;

  ierr = PetscNewLog(tao,&tron);CHKERRQ(ierr);
  tao->data = (void*)tron;

  /* Override default settings (unless already changed) */
  if (!tao->max_it_changed) tao->max_it = 50;

#if defined(PETSC_USE_REAL_SINGLE)
  if (!tao->steptol_changed) tao->steptol = 1.0e-6;
#else
  if (!tao->steptol_changed) tao->steptol = 1.0e-12;
#endif

  if (!tao->trust0_changed) tao->trust0 = 1.0;

  /* Initialize pointers and variables */
  tron->n            = 0;
  tron->maxgpits     = 3;
  tron->pg_ftol      = 0.001;

  tron->eta1         = 1.0e-4;
  tron->eta2         = 0.25;
  tron->eta3         = 0.50;
  tron->eta4         = 0.90;

  tron->sigma1       = 0.5;
  tron->sigma2       = 2.0;
  tron->sigma3       = 4.0;

  tron->gp_iterates  = 0; /* Cumulative number */
  tron->total_gp_its = 0;
  tron->n_free       = 0;

  tron->DXFree=NULL;
  tron->R=NULL;
  tron->X_New=NULL;
  tron->G_New=NULL;
  tron->Work=NULL;
  tron->Free_Local=NULL;
  tron->H_sub=NULL;
  tron->Hpre_sub=NULL;
  tao->subset_type = TAO_SUBSET_SUBVEC;

  ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
  ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr);
  ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
  ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);

  ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr);
  ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr);
  ierr = KSPSetType(tao->ksp,KSPSTCG);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}
开发者ID:plguhur,项目名称:petsc,代码行数:66,代码来源:tron.c

示例10: main

int main(int argc,char **args)
{
  KSP            subksp;
  Mat            A,subA;
  Vec            x,b,u,subb,subx,subu;
  PetscViewer    fd;
  char           file[PETSC_MAX_PATH_LEN];
  PetscBool      flg;
  PetscErrorCode ierr;
  PetscInt       i,m,n,its;
  PetscReal      norm;
  PetscMPIInt    rank,size;
  MPI_Comm       comm,subcomm;
  PetscSubcomm   psubcomm;
  PetscInt       nsubcomm=1,id;
  PetscScalar    *barray,*xarray,*uarray,*array,one=1.0;
  PetscInt       type=1;

  PetscInitialize(&argc,&args,(char*)0,help);
  /* Load the matrix */
  ierr = PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
  if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr);

  /* Load the matrix; then destroy the viewer.*/
  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatLoad(A,fd);CHKERRQ(ierr);
  ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr);

  ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);

  /* Create rhs vector b */
  ierr = MatGetLocalSize(A,&m,NULL);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecSetSizes(b,m,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);
  ierr = VecSet(b,one);CHKERRQ(ierr);

  ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&u);CHKERRQ(ierr);
  ierr = VecSet(x,0.0);CHKERRQ(ierr);

  /* Test MatGetMultiProcBlock() */
  ierr = PetscOptionsGetInt(NULL,"-nsubcomm",&nsubcomm,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsGetInt(NULL,"-subcomm_type",&type,NULL);CHKERRQ(ierr);

  ierr = PetscSubcommCreate(comm,&psubcomm);CHKERRQ(ierr);
  ierr = PetscSubcommSetNumber(psubcomm,nsubcomm);CHKERRQ(ierr);
  if (type == PETSC_SUBCOMM_GENERAL) { /* user provides color, subrank and duprank */
    PetscMPIInt color,subrank,duprank,subsize;
    duprank = size-1 - rank;
    subsize = size/nsubcomm;
    if (subsize*nsubcomm != size) SETERRQ2(comm,PETSC_ERR_SUP,"This example requires nsubcomm %D divides nproc %D",nsubcomm,size);
    color   = duprank/subsize;
    subrank = duprank - color*subsize;
    ierr    = PetscSubcommSetTypeGeneral(psubcomm,color,subrank,duprank);CHKERRQ(ierr);
  } else if (type == PETSC_SUBCOMM_CONTIGUOUS) {
    ierr = PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr);
  } else if (type == PETSC_SUBCOMM_INTERLACED) {
    ierr = PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_INTERLACED);CHKERRQ(ierr);
  } else SETERRQ1(psubcomm->parent,PETSC_ERR_SUP,"PetscSubcommType %D is not supported yet",type);
  subcomm = psubcomm->comm;

  ierr = PetscOptionsHasName(NULL, "-subcomm_view", &flg);CHKERRQ(ierr);
  if (flg) {
    PetscMPIInt subsize,subrank,duprank;
    ierr = MPI_Comm_size((MPI_Comm)subcomm,&subsize);CHKERRQ(ierr);
    ierr = MPI_Comm_rank((MPI_Comm)subcomm,&subrank);CHKERRQ(ierr);
    ierr = MPI_Comm_rank((MPI_Comm)psubcomm->dupparent,&duprank);CHKERRQ(ierr);

    ierr = PetscSynchronizedPrintf(comm,"[%D], color %D, sub-size %D, sub-rank %D, duprank %D\n",rank,psubcomm->color,subsize,subrank,duprank);
    ierr = PetscSynchronizedFlush(comm);CHKERRQ(ierr);
  }

  /* Create subA */
  ierr = MatGetMultiProcBlock(A,subcomm,MAT_INITIAL_MATRIX,&subA);CHKERRQ(ierr);

  /* Create sub vectors without arrays. Place b's and x's local arrays into subb and subx */
  ierr = MatGetLocalSize(subA,&m,&n);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&subb);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&subx);CHKERRQ(ierr);
  ierr = VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&subu);CHKERRQ(ierr);

  ierr = VecGetArray(b,&barray);CHKERRQ(ierr);
  ierr = VecGetArray(x,&xarray);CHKERRQ(ierr);
  ierr = VecGetArray(u,&uarray);CHKERRQ(ierr);
  ierr = VecPlaceArray(subb,barray);CHKERRQ(ierr);
  ierr = VecPlaceArray(subx,xarray);CHKERRQ(ierr);
  ierr = VecPlaceArray(subu,uarray);CHKERRQ(ierr);

  /* Create linear solvers associated with subA */
  ierr = KSPCreate(subcomm,&subksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(subksp,subA,subA,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(subksp);CHKERRQ(ierr);

  /* Solve sub systems */
  ierr = KSPSolve(subksp,subb,subx);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(subksp,&its);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:feelpp,项目名称:debian-petsc,代码行数:101,代码来源:ex37.c

示例11: BSSCR_DRIVER_flex


//.........这里部分代码省略.........
    useNormInfStoppingConditions = PETSC_FALSE;
    PetscOptionsGetTruth( PETSC_NULL ,"-A11_use_norm_inf_stopping_condition", &useNormInfStoppingConditions, &found );
    if(useNormInfStoppingConditions) 
        BSSCR_KSPSetNormInfConvergenceTest( ksp_inner ); 

    useNormInfMonitor = PETSC_FALSE; 
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor) 
        KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfMonitor, PETSC_NULL, PETSC_NULL );
        
    useNormInfMonitor = PETSC_FALSE; 
    PetscOptionsGetTruth( PETSC_NULL, "-A11_ksp_norm_inf_to_norm_2_monitor", &useNormInfMonitor, &found );
    if(useNormInfMonitor) 
        KSPMonitorSet( ksp_inner, BSSCR_KSPNormInfToNorm2Monitor, PETSC_NULL, PETSC_NULL );

    /* create right hand side */
    /* h_hat = G'*inv(K)*f - h */
    MatGetVecs(K,PETSC_NULL,&t);
    MatGetVecs( S, PETSC_NULL, &h_hat );

    KSPSetOptionsPrefix( ksp_inner, "A11_" );
    KSPSetFromOptions( ksp_inner );

    KSPSolve(ksp_inner,f,t);/* t=f/K */
    //bsscr_writeVec( t, "ts", "Writing t vector");
    MatMult(D,t,h_hat);/* G'*t */
    VecAXPY(h_hat, -1.0, h);/* h_hat = h_hat - h */
    Stg_VecDestroy(&t);
    //bsscr_writeVec( h_hat, "h_hat", "Writing h_hat Vector in Solver");
    //MatSchurApplyReductionToVecFromBlock( S, stokes_b, h_hat );
        
    /* create solver for S p = h_hat */

    KSPCreate( PETSC_COMM_WORLD, &ksp_S );
    KSPSetOptionsPrefix( ksp_S, "scr_");
    Stg_KSPSetOperators( ksp_S, S,S, SAME_NONZERO_PATTERN );
    KSPSetType( ksp_S, "cg" );
        
    /* Build preconditioner for S */
    KSPGetPC( ksp_S, &pc_S );
    BSSCR_BSSCR_StokesCreatePCSchur2( K,G,D,C,approxS,pc_S,sym, bsscrp_self );

    KSPSetFromOptions(ksp_S);

    /* Set specific monitor test */
    KSPGetTolerances( ksp_S, PETSC_NULL, PETSC_NULL, PETSC_NULL, &max_it );
    //BSSCR_KSPLogSetMonitor( ksp_S, max_it, &monitor_index );
        
    /* Pressure / Velocity Solve */   
    scrSolveTime = MPI_Wtime();
    PetscPrintf( PETSC_COMM_WORLD, "\t* Pressure / Velocity Solve \n");

    usePreviousGuess = PETSC_FALSE; 
    if(been_here)
        PetscOptionsGetTruth( PETSC_NULL, "-scr_use_previous_guess", &usePreviousGuess, &found );
    
    if(usePreviousGuess) {   /* Note this should actually look at checkpoint information */
        KSPSetInitialGuessNonzero( ksp_S, PETSC_TRUE );
    }
    else {
        KSPSetInitialGuessNonzero( ksp_S, PETSC_FALSE );
    }
    
    //KSPSetRelativeRhsConvergenceTest( ksp_S );

    useNormInfStoppingConditions = PETSC_FALSE;
开发者ID:OlympusMonds,项目名称:EarthByte_Underworld,代码行数:67,代码来源:bsscr-driver.c

示例12: main

int main(int argc, char** argv)
{
    DM da;
    PetscErrorCode ierr;
    Vec x, rhs;
    Mat A, jac;
    ierr = PetscInitialize(&argc, &argv, NULL, NULL);
    CHKERRQ(ierr);
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Laplacian in 2D", "");
    CHKERRQ(ierr);
    ierr = PetscOptionsEnd();
    CHKERRQ(ierr);
    ierr = HpddmRegisterKSP();
    CHKERRQ(ierr);
    MPI_Barrier(PETSC_COMM_WORLD);
    double time = MPI_Wtime();
    ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 10, 10, PETSC_DECIDE, PETSC_DECIDE, 1, 1,
                        0, 0, &da);
    CHKERRQ(ierr);
    ierr = DMSetFromOptions(da);
    CHKERRQ(ierr);
    ierr = DMSetUp(da);
    CHKERRQ(ierr);
    ierr = DMCreateGlobalVector(da, &rhs);
    CHKERRQ(ierr);
    ierr = DMCreateGlobalVector(da, &x);
    CHKERRQ(ierr);
    ierr = DMCreateMatrix(da, &A);
    CHKERRQ(ierr);
    ierr = DMCreateMatrix(da, &jac);
    CHKERRQ(ierr);
    ierr = ComputeMatrix(da, jac, A);
    CHKERRQ(ierr);
    MPI_Barrier(PETSC_COMM_WORLD);
    time = MPI_Wtime() - time;
    ierr = PetscPrintf(PETSC_COMM_WORLD, "--- Mat assembly = %f\n", time);
    CHKERRQ(ierr);
    MPI_Barrier(PETSC_COMM_WORLD);
    time = MPI_Wtime();
    KSP ksp;
    ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);
    CHKERRQ(ierr);
    ierr = KSPSetDM(ksp, da);
    CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);
    CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp, A, A);
    CHKERRQ(ierr);
    ierr = KSPSetDMActive(ksp, PETSC_FALSE);
    CHKERRQ(ierr);
    ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
    CHKERRQ(ierr);
    ierr = KSPSetUp(ksp);
    CHKERRQ(ierr);
    MPI_Barrier(PETSC_COMM_WORLD);
    time = MPI_Wtime() - time;
    ierr = PetscPrintf(PETSC_COMM_WORLD, "--- PC setup = %f\n", time);
    CHKERRQ(ierr);
    PetscScalar nus[SIZE_ARRAY_NU] = {0.1, 10.0, 0.001, 100.0};
    float t_time[SIZE_ARRAY_NU];
    int t_its[SIZE_ARRAY_NU];
    int i, j;
    for (j = 0; j < 2; ++j) {
        {
            if (j == 1) {
                ierr = KSPSetType(ksp, "hpddm");
                CHKERRQ(ierr);
                ierr = KSPSetFromOptions(ksp);
                CHKERRQ(ierr);
                ierr = VecZeroEntries(x);
                CHKERRQ(ierr);
            }
            ierr = KSPSolve(ksp, rhs, x);
            CHKERRQ(ierr);
            if (j == 1) {
                const HpddmOption* const opt = HpddmOptionGet();
                int previous = HpddmOptionVal(opt, "krylov_method");
                if (previous == HPDDM_KRYLOV_METHOD_GCRODR || previous == HPDDM_KRYLOV_METHOD_BGCRODR) HpddmDestroyRecycling();
            }
        }
        for (i = 0; i < SIZE_ARRAY_NU; ++i) {
            ierr = VecZeroEntries(x);
            CHKERRQ(ierr);
            ierr = ComputeRHS(da, rhs, nus[i]);
            CHKERRQ(ierr);
            MPI_Barrier(PETSC_COMM_WORLD);
            time = MPI_Wtime();
            ierr = KSPSolve(ksp, rhs, x);
            CHKERRQ(ierr);
            MPI_Barrier(PETSC_COMM_WORLD);
            t_time[i] = MPI_Wtime() - time;
            PetscInt its;
            ierr = KSPGetIterationNumber(ksp, &its);
            CHKERRQ(ierr);
            t_its[i] = its;
            ierr = ComputeError(A, rhs, x);
            CHKERRQ(ierr);
        }
        for (i = 0; i < SIZE_ARRAY_NU; ++i) {
            ierr = PetscPrintf(PETSC_COMM_WORLD, "%d\t%d\t%f\n", i + 1, t_its[i], t_time[i]);
//.........这里部分代码省略.........
开发者ID:hpddm,项目名称:hpddm,代码行数:101,代码来源:ex32.c

示例13: main

int main(int argc,char **args)
{
  Mat            C;
  PetscMPIInt    rank,size;
  PetscInt       i,m = 5,N,start,end,M,its;
  PetscScalar    val,Ke[16],r[4];
  PetscReal      x,y,h,norm,tol=1.e-14;
  PetscErrorCode ierr;
  PetscInt       idx[4],count,*rows;
  Vec            u,ustar,b;
  KSP            ksp;

  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);
  ierr  = MatSetUp(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 */
    x = h*(i % m); y = h*(i/m);
    /* 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 = VecDuplicate(b,&ustar);CHKERRQ(ierr);
  ierr = VecSet(u,0.0);CHKERRQ(ierr);
  ierr = VecSet(b,0.0);CHKERRQ(ierr);

  /* Assemble right-hand-side vector */
  for (i=start; i<end; i++) {
    /* location of lower left corner of element */
    x = h*(i % m); y = h*(i/m);
    /* 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   = FormElementRhs(x,y,h*h,r);CHKERRQ(ierr);
    ierr   = VecSetValues(b,4,idx,r,ADD_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

  /* Modify matrix and right-hand-side for Dirichlet boundary conditions */
  ierr = PetscMalloc1(4*m,&rows);CHKERRQ(ierr);
  for (i=0; i<m+1; i++) {
    rows[i]          = i; /* bottom */
    rows[3*m - 1 +i] = m*(m+1) + i; /* top */
  }
  count = m+1; /* left side */
  for (i=m+1; i<m*(m+1); i+= m+1) rows[count++] = i;

  count = 2*m; /* left side */
  for (i=2*m+1; i<m*(m+1); i+= m+1) rows[count++] = i;
  for (i=0; i<4*m; i++) {
    x    = h*(rows[i] % (m+1)); y = h*(rows[i]/(m+1));
    val  = y;
    ierr = VecSetValues(u,1,&rows[i],&val,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(b,1,&rows[i],&val,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatZeroRows(C,4*m,rows,1.0,0,0);CHKERRQ(ierr);

  ierr = PetscFree(rows);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

  { Mat A;
    ierr = MatConvert(C,MATSAME,MAT_INITIAL_MATRIX,&A);CHKERRQ(ierr);
    ierr = MatDestroy(&C);CHKERRQ(ierr);
    ierr = MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);CHKERRQ(ierr);
    ierr = MatDestroy(&A);CHKERRQ(ierr);
  }

  /* Solve linear system */
  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,C,C,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:ZJLi2013,项目名称:petsc,代码行数:101,代码来源:ex3.c

示例14: main

int main(int argc,char **args)
{
  Mat            C;
  PetscErrorCode ierr;
  PetscInt       i,m = 2,N,M,its,idx[4],count,*rows;
  PetscScalar    val,Ke[16],r[4];
  PetscReal      x,y,h,norm,tol=1.e-14;
  Vec            u,ustar,b;
  KSP            ksp;

  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 */

  /* create stiffness matrix */
  ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,9,NULL,&C);CHKERRQ(ierr);
  ierr = MatSetUp(C);CHKERRQ(ierr);

  /* forms the element stiffness for the Laplacian */
  ierr = FormElementStiffness(h*h,Ke);CHKERRQ(ierr);
  for (i=0; i<M; i++) {
    /* location of lower left corner of element */
    x = h*(i % m); y = h*(i/m);
    /* 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 */

  ierr = VecCreateSeq(PETSC_COMM_SELF,N,&u);CHKERRQ(ierr);
  ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
  ierr = VecDuplicate(b,&ustar);CHKERRQ(ierr);
  ierr = VecSet(u,0.0);CHKERRQ(ierr);
  ierr = VecSet(b,0.0);CHKERRQ(ierr);

  for (i=0; i<M; i++) {
    /* location of lower left corner of element */
    x = h*(i % m); y = h*(i/m);
    /* 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   = FormElementRhs(x,y,h*h,r);CHKERRQ(ierr);
    ierr   = VecSetValues(b,4,idx,r,ADD_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

  /* modify matrix and rhs for Dirichlet boundary conditions */
  ierr = PetscMalloc1((4*m+1),&rows);CHKERRQ(ierr);
  for (i=0; i<m+1; i++) {
    rows[i]          = i; /* bottom */
    rows[3*m - 1 +i] = m*(m+1) + i; /* top */
  }
  count = m+1; /* left side */
  for (i=m+1; i<m*(m+1); i+= m+1) rows[count++] = i;

  count = 2*m; /* left side */
  for (i=2*m+1; i<m*(m+1); i+= m+1) rows[count++] = i;

  for (i=0; i<4*m; i++) {
    x    = h*(rows[i] % (m+1)); y = h*(rows[i]/(m+1));
    val  = y;
    ierr = VecSetValues(u,1,&rows[i],&val,INSERT_VALUES);CHKERRQ(ierr);
    ierr = VecSetValues(b,1,&rows[i],&val,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = MatZeroRows(C,4*m,rows,1.0,0,0);CHKERRQ(ierr);

  ierr = PetscFree(rows);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(u);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(u);CHKERRQ(ierr);
  ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(b);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);
  ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);

  /* check error */
  for (i=0; i<N; i++) {
    x    = h*(i % (m+1)); y = h*(i/(m+1));
    val  = y;
    ierr = VecSetValues(ustar,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr);
  }
  ierr = VecAssemblyBegin(ustar);CHKERRQ(ierr);
  ierr = VecAssemblyEnd(ustar);CHKERRQ(ierr);

  ierr = VecAXPY(u,-1.0,ustar);CHKERRQ(ierr);
  ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr);
  ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
  if (norm > tol) {
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G Iterations %D\n",norm*h,its);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:ZJLi2013,项目名称:petsc,代码行数:101,代码来源:ex4.c

示例15: PCSetUp_Redundant

static PetscErrorCode PCSetUp_Redundant(PC pc)
{
  PC_Redundant   *red = (PC_Redundant*)pc->data;
  PetscErrorCode ierr;
  PetscInt       mstart,mend,mlocal,M;
  PetscMPIInt    size;
  MPI_Comm       comm,subcomm;
  Vec            x;
  const char     *prefix;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
  
  /* if pmatrix set by user is sequential then we do not need to gather the parallel matrix */
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  if (size == 1) red->useparallelmat = PETSC_FALSE;

  if (!pc->setupcalled) {
    PetscInt mloc_sub;
    if (!red->psubcomm) {
      ierr = PetscSubcommCreate(comm,&red->psubcomm);CHKERRQ(ierr);
      ierr = PetscSubcommSetNumber(red->psubcomm,red->nsubcomm);CHKERRQ(ierr);
      ierr = PetscSubcommSetType(red->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);CHKERRQ(ierr);
      /* enable runtime switch of psubcomm type, e.g., '-psubcomm_type interlaced */
      ierr = PetscSubcommSetFromOptions(red->psubcomm);CHKERRQ(ierr);
      ierr = PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));CHKERRQ(ierr);

      /* create a new PC that processors in each subcomm have copy of */
      subcomm = PetscSubcommChild(red->psubcomm);

      ierr = KSPCreate(subcomm,&red->ksp);CHKERRQ(ierr);
      ierr = KSPSetErrorIfNotConverged(red->ksp,pc->erroriffailure);CHKERRQ(ierr);
      ierr = PetscObjectIncrementTabLevel((PetscObject)red->ksp,(PetscObject)pc,1);CHKERRQ(ierr);
      ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)red->ksp);CHKERRQ(ierr);
      ierr = KSPSetType(red->ksp,KSPPREONLY);CHKERRQ(ierr);
      ierr = KSPGetPC(red->ksp,&red->pc);CHKERRQ(ierr);
      ierr = PCSetType(red->pc,PCLU);CHKERRQ(ierr);

      ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr);
      ierr = KSPSetOptionsPrefix(red->ksp,prefix);CHKERRQ(ierr);
      ierr = KSPAppendOptionsPrefix(red->ksp,"redundant_");CHKERRQ(ierr);
    } else {
      subcomm = PetscSubcommChild(red->psubcomm);
    }

    if (red->useparallelmat) {
      /* grab the parallel matrix and put it into processors of a subcomminicator */
      ierr = MatCreateRedundantMatrix(pc->pmat,red->psubcomm->n,subcomm,MAT_INITIAL_MATRIX,&red->pmats);CHKERRQ(ierr);
      ierr = KSPSetOperators(red->ksp,red->pmats,red->pmats);CHKERRQ(ierr);
       
      /* get working vectors xsub and ysub */
      ierr = MatCreateVecs(red->pmats,&red->xsub,&red->ysub);CHKERRQ(ierr);

      /* create working vectors xdup and ydup.
       xdup concatenates all xsub's contigously to form a mpi vector over dupcomm  (see PetscSubcommCreate_interlaced())
       ydup concatenates all ysub and has empty local arrays because ysub's arrays will be place into it.
       Note: we use communicator dupcomm, not PetscObjectComm((PetscObject)pc)! */
      ierr = MatGetLocalSize(red->pmats,&mloc_sub,NULL);CHKERRQ(ierr);
      ierr = VecCreateMPI(PetscSubcommContiguousParent(red->psubcomm),mloc_sub,PETSC_DECIDE,&red->xdup);CHKERRQ(ierr);
      ierr = VecCreateMPIWithArray(PetscSubcommContiguousParent(red->psubcomm),1,mloc_sub,PETSC_DECIDE,NULL,&red->ydup);CHKERRQ(ierr);

      /* create vecscatters */
      if (!red->scatterin) { /* efficiency of scatterin is independent from psubcomm_type! */
        IS       is1,is2;
        PetscInt *idx1,*idx2,i,j,k;

        ierr = MatCreateVecs(pc->pmat,&x,0);CHKERRQ(ierr);
        ierr = VecGetSize(x,&M);CHKERRQ(ierr);
        ierr = VecGetOwnershipRange(x,&mstart,&mend);CHKERRQ(ierr);
        mlocal = mend - mstart;
        ierr = PetscMalloc2(red->psubcomm->n*mlocal,&idx1,red->psubcomm->n*mlocal,&idx2);CHKERRQ(ierr);
        j    = 0;
        for (k=0; k<red->psubcomm->n; k++) {
          for (i=mstart; i<mend; i++) {
            idx1[j]   = i;
            idx2[j++] = i + M*k;
          }
        }
        ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx1,PETSC_COPY_VALUES,&is1);CHKERRQ(ierr);
        ierr = ISCreateGeneral(comm,red->psubcomm->n*mlocal,idx2,PETSC_COPY_VALUES,&is2);CHKERRQ(ierr);
        ierr = VecScatterCreate(x,is1,red->xdup,is2,&red->scatterin);CHKERRQ(ierr);
        ierr = ISDestroy(&is1);CHKERRQ(ierr);
        ierr = ISDestroy(&is2);CHKERRQ(ierr);

        /* Impl below is good for PETSC_SUBCOMM_INTERLACED (no inter-process communication) and PETSC_SUBCOMM_CONTIGUOUS (communication within subcomm) */
        ierr = ISCreateStride(comm,mlocal,mstart+ red->psubcomm->color*M,1,&is1);CHKERRQ(ierr);
        ierr = ISCreateStride(comm,mlocal,mstart,1,&is2);CHKERRQ(ierr);
        ierr = VecScatterCreate(red->xdup,is1,x,is2,&red->scatterout);CHKERRQ(ierr);
        ierr = ISDestroy(&is1);CHKERRQ(ierr);
        ierr = ISDestroy(&is2);CHKERRQ(ierr);
        ierr = PetscFree2(idx1,idx2);CHKERRQ(ierr);
        ierr = VecDestroy(&x);CHKERRQ(ierr);
      }
    } else { /* !red->useparallelmat */
      ierr = KSPSetOperators(red->ksp,pc->mat,pc->pmat);CHKERRQ(ierr);
    }
  } else { /* pc->setupcalled */
    if (red->useparallelmat) {
      MatReuse       reuse;
      /* grab the parallel matrix and put it into processors of a subcomminicator */
//.........这里部分代码省略.........
开发者ID:plguhur,项目名称:petsc,代码行数:101,代码来源:redundant.c


注:本文中的KSPCreate函数示例由纯净天空整理自Github/MSDocs等开源代码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。