• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

C++ PetscMalloc1函数代码示例

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

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



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

示例1: DMPlexComputeAnchorAdjacencies

/* get adjacencies due to point-to-point constraints that can't be found with DMPlexGetAdjacency() */
static PetscErrorCode DMPlexComputeAnchorAdjacencies(DM dm, PetscSection section, PetscSection sectionGlobal, PetscSection *anchorSectionAdj, PetscInt *anchorAdj[])
{
  PetscInt       pStart, pEnd;
  PetscSection   adjSec, aSec;
  IS             aIS;
  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = PetscSectionCreate(PetscObjectComm((PetscObject)section),&adjSec);CHKERRQ(ierr);
  ierr = PetscSectionGetChart(section,&pStart,&pEnd);CHKERRQ(ierr);
  ierr = PetscSectionSetChart(adjSec,pStart,pEnd);CHKERRQ(ierr);

  ierr = DMPlexGetAnchors(dm,&aSec,&aIS);CHKERRQ(ierr);
  if (aSec) {
    const PetscInt *anchors;
    PetscInt       p, q, a, aSize, *offsets, aStart, aEnd, *inverse, iSize, *adj, adjSize;
    PetscInt       *tmpAdjP = NULL, *tmpAdjQ = NULL;
    PetscSection   inverseSec;

    /* invert the constraint-to-anchor map */
    ierr = PetscSectionCreate(PetscObjectComm((PetscObject)aSec),&inverseSec);CHKERRQ(ierr);
    ierr = PetscSectionSetChart(inverseSec,pStart,pEnd);CHKERRQ(ierr);
    ierr = ISGetLocalSize(aIS, &aSize);CHKERRQ(ierr);
    ierr = ISGetIndices(aIS, &anchors);CHKERRQ(ierr);

    for (p = 0; p < aSize; p++) {
      PetscInt a = anchors[p];

      ierr = PetscSectionAddDof(inverseSec,a,1);CHKERRQ(ierr);
    }
    ierr = PetscSectionSetUp(inverseSec);CHKERRQ(ierr);
    ierr = PetscSectionGetStorageSize(inverseSec,&iSize);CHKERRQ(ierr);
    ierr = PetscMalloc1(iSize,&inverse);CHKERRQ(ierr);
    ierr = PetscCalloc1(pEnd-pStart,&offsets);CHKERRQ(ierr);
    ierr = PetscSectionGetChart(aSec,&aStart,&aEnd);CHKERRQ(ierr);
    for (p = aStart; p < aEnd; p++) {
      PetscInt dof, off;

      ierr = PetscSectionGetDof(aSec, p, &dof);CHKERRQ(ierr);
      ierr = PetscSectionGetOffset(aSec, p, &off);CHKERRQ(ierr);

      for (q = 0; q < dof; q++) {
        PetscInt iOff;

        a = anchors[off + q];
        ierr = PetscSectionGetOffset(inverseSec, a, &iOff);CHKERRQ(ierr);
        inverse[iOff + offsets[a-pStart]++] = p;
      }
    }
    ierr = ISRestoreIndices(aIS, &anchors);CHKERRQ(ierr);
    ierr = PetscFree(offsets);CHKERRQ(ierr);

    /* construct anchorAdj and adjSec
     *
     * loop over anchors:
     *   construct anchor adjacency
     *   loop over constrained:
     *     construct constrained adjacency
     *     if not in anchor adjacency, add to dofs
     * setup adjSec, allocate anchorAdj
     * loop over anchors:
     *   construct anchor adjacency
     *   loop over constrained:
     *     construct constrained adjacency
     *     if not in anchor adjacency
     *       if not already in list, put in list
     *   sort, unique, reduce dof count
     * optional: compactify
     */
    for (p = pStart; p < pEnd; p++) {
      PetscInt iDof, iOff, i, r, s, numAdjP = PETSC_DETERMINE;

      ierr = PetscSectionGetDof(inverseSec,p,&iDof);CHKERRQ(ierr);
      if (!iDof) continue;
      ierr = PetscSectionGetOffset(inverseSec,p,&iOff);CHKERRQ(ierr);
      ierr = DMPlexGetAdjacency(dm,p,&numAdjP,&tmpAdjP);CHKERRQ(ierr);
      for (i = 0; i < iDof; i++) {
        PetscInt iNew = 0, qAdj, qAdjDof, qAdjCDof, numAdjQ = PETSC_DETERMINE;

        q = inverse[iOff + i];
        ierr = DMPlexGetAdjacency(dm,q,&numAdjQ,&tmpAdjQ);CHKERRQ(ierr);
        for (r = 0; r < numAdjQ; r++) {
          qAdj = tmpAdjQ[r];
          if ((qAdj < pStart) || (qAdj >= pEnd)) continue;
          for (s = 0; s < numAdjP; s++) {
            if (qAdj == tmpAdjP[s]) break;
          }
          if (s < numAdjP) continue;
          ierr  = PetscSectionGetDof(section,qAdj,&qAdjDof);CHKERRQ(ierr);
          ierr  = PetscSectionGetConstraintDof(section,qAdj,&qAdjCDof);CHKERRQ(ierr);
          iNew += qAdjDof - qAdjCDof;
        }
        ierr = PetscSectionAddDof(adjSec,p,iNew);CHKERRQ(ierr);
      }
    }

    ierr = PetscSectionSetUp(adjSec);CHKERRQ(ierr);
    ierr = PetscSectionGetStorageSize(adjSec,&adjSize);CHKERRQ(ierr);
//.........这里部分代码省略.........
开发者ID:jgoldfar,项目名称:petsc,代码行数:101,代码来源:plexpreallocate.c


示例2: MatSetUpMultiply_MPIAIJ

PetscErrorCode MatSetUpMultiply_MPIAIJ(Mat mat)
{
  Mat_MPIAIJ     *aij = (Mat_MPIAIJ*)mat->data;
  Mat_SeqAIJ     *B   = (Mat_SeqAIJ*)(aij->B->data);
  PetscErrorCode ierr;
  PetscInt       i,j,*aj = B->j,ec = 0,*garray;
  IS             from,to;
  Vec            gvec;
#if defined(PETSC_USE_CTABLE)
  PetscTable         gid1_lid1;
  PetscTablePosition tpos;
  PetscInt           gid,lid;
#else
  PetscInt N = mat->cmap->N,*indices;
#endif

  PetscFunctionBegin;
  if (!aij->garray) {
#if defined(PETSC_USE_CTABLE)
    /* use a table */
    ierr = PetscTableCreate(aij->B->rmap->n,mat->cmap->N+1,&gid1_lid1);CHKERRQ(ierr);
    for (i=0; i<aij->B->rmap->n; i++) {
      for (j=0; j<B->ilen[i]; j++) {
        PetscInt data,gid1 = aj[B->i[i] + j] + 1;
        ierr = PetscTableFind(gid1_lid1,gid1,&data);CHKERRQ(ierr);
        if (!data) {
          /* one based table */
          ierr = PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES);CHKERRQ(ierr);
        }
      }
    }
    /* form array of columns we need */
    ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr);
    ierr = PetscTableGetHeadPosition(gid1_lid1,&tpos);CHKERRQ(ierr);
    while (tpos) {
      ierr = PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid);CHKERRQ(ierr);
      gid--;
      lid--;
      garray[lid] = gid;
    }
    ierr = PetscSortInt(ec,garray);CHKERRQ(ierr); /* sort, and rebuild */
    ierr = PetscTableRemoveAll(gid1_lid1);CHKERRQ(ierr);
    for (i=0; i<ec; i++) {
      ierr = PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES);CHKERRQ(ierr);
    }
    /* compact out the extra columns in B */
    for (i=0; i<aij->B->rmap->n; i++) {
      for (j=0; j<B->ilen[i]; j++) {
        PetscInt gid1 = aj[B->i[i] + j] + 1;
        ierr = PetscTableFind(gid1_lid1,gid1,&lid);CHKERRQ(ierr);
        lid--;
        aj[B->i[i] + j] = lid;
      }
    }
    aij->B->cmap->n = aij->B->cmap->N = ec;
    aij->B->cmap->bs = 1;

    ierr = PetscLayoutSetUp((aij->B->cmap));CHKERRQ(ierr);
    ierr = PetscTableDestroy(&gid1_lid1);CHKERRQ(ierr);
#else
    /* Make an array as long as the number of columns */
    /* mark those columns that are in aij->B */
    ierr = PetscCalloc1(N+1,&indices);CHKERRQ(ierr);
    for (i=0; i<aij->B->rmap->n; i++) {
      for (j=0; j<B->ilen[i]; j++) {
        if (!indices[aj[B->i[i] + j]]) ec++;
        indices[aj[B->i[i] + j]] = 1;
      }
    }

    /* form array of columns we need */
    ierr = PetscMalloc1(ec+1,&garray);CHKERRQ(ierr);
    ec   = 0;
    for (i=0; i<N; i++) {
      if (indices[i]) garray[ec++] = i;
    }

    /* make indices now point into garray */
    for (i=0; i<ec; i++) {
      indices[garray[i]] = i;
    }

    /* compact out the extra columns in B */
    for (i=0; i<aij->B->rmap->n; i++) {
      for (j=0; j<B->ilen[i]; j++) {
        aj[B->i[i] + j] = indices[aj[B->i[i] + j]];
      }
    }
    aij->B->cmap->n = aij->B->cmap->N = ec;
    aij->B->cmap->bs = 1;

    ierr = PetscLayoutSetUp((aij->B->cmap));CHKERRQ(ierr);
    ierr = PetscFree(indices);CHKERRQ(ierr);
#endif
  } else {
    garray = aij->garray;
  }

  if (!aij->lvec) {
    /* create local vector that is used to scatter into */
//.........这里部分代码省略.........
开发者ID:petsc,项目名称:petsc,代码行数:101,代码来源:mmaij.c


示例3: main


//.........这里部分代码省略.........
  VecAssemblyEnd(fin);
/*    VecView(fin,PETSC_VIEWER_STDOUT_WORLD); */


  VecGetArray(fin,&x_arr);
  VecGetArray(fout1,&z_arr);
  VecGetArray(fout,&y_arr);

  fplan=fftw_mpi_plan_dft_r2c_3d(N0,N1,N2,(double*)x_arr,(fftw_complex*)y_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);
  bplan=fftw_mpi_plan_dft_c2r_3d(N0,N1,N2,(fftw_complex*)y_arr,(double*)z_arr,PETSC_COMM_WORLD,FFTW_ESTIMATE);

  fftw_execute(fplan);
  fftw_execute(bplan);

  VecRestoreArray(fin,&x_arr);
  VecRestoreArray(fout1,&z_arr);
  VecRestoreArray(fout,&y_arr);


/*    a = 1.0/(PetscReal)N_factor; */
/*    ierr = VecScale(fout1,a);CHKERRQ(ierr); */
  VecCreate(PETSC_COMM_WORLD,&ini);
  VecCreate(PETSC_COMM_WORLD,&final);
  VecSetSizes(ini,local_n0*N1*N2,N_factor);
  VecSetSizes(final,local_n0*N1*N2,N_factor);
/*    VecSetSizes(ini,PETSC_DECIDE,N_factor); */
/*    VecSetSizes(final,PETSC_DECIDE,N_factor); */
  VecSetFromOptions(ini);
  VecSetFromOptions(final);

  if (N2%2==0) NM=N2+2;
  else NM=N2+1;

  ierr = VecGetOwnershipRange(fin,&low,NULL);
  printf("The local index is %d from %d\n",low,rank);
  ierr = PetscMalloc1(local_n0*N1*N2,&indx3);
  ierr = PetscMalloc1(local_n0*N1*N2,&indx4);
  for (i=0; i<local_n0; i++) {
    for (j=0;j<N1;j++) {
      for (k=0;k<N2;k++) {
        tempindx  = i*N1*N2 + j*N2 + k;
        tempindx1 = i*N1*NM + j*NM + k;

        indx3[tempindx]=local_0_start*N1*N2+tempindx;
        indx4[tempindx]=low+tempindx1;
      }
      /*          printf("index3 %d from proc %d is \n",indx3[tempindx],rank); */
      /*          printf("index4 %d from proc %d is \n",indx4[tempindx],rank); */
    }
  }
  VecGetValues(fin,local_n0*N1*N2,indx4,x_arr);
  VecSetValues(ini,local_n0*N1*N2,indx3,x_arr,INSERT_VALUES);
  VecAssemblyBegin(ini);
  VecAssemblyEnd(ini);

  VecGetValues(fout1,local_n0*N1*N2,indx4,y_arr);
  VecSetValues(final,local_n0*N1*N2,indx3,y_arr,INSERT_VALUES);
  VecAssemblyBegin(final);
  VecAssemblyEnd(final);

  printf("The local index value is %ld from %d",local_n0*N1*N2,rank);
/*
  for (i=0;i<N0;i++) {
     for (j=0;j<N1;j++) {
        indx=i*N1*NM+j*NM;
        ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx1);
        indx=i*N1*N2+j*N2;
        ISCreateStride(PETSC_COMM_WORLD,N2,indx,1,&indx2);
        VecScatterCreate(fin,indx1,ini,indx2,&vecscat);
        VecScatterBegin(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
        VecScatterEnd(vecscat,fin,ini,INSERT_VALUES,SCATTER_FORWARD);
        VecScatterCreate(fout1,indx1,final,indx2,&vecscat1);
        VecScatterBegin(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
        VecScatterEnd(vecscat1,fout1,final,INSERT_VALUES,SCATTER_FORWARD);
     }
  }
*/
  a    = 1.0/(PetscReal)N_factor;
  ierr = VecScale(fout1,a);CHKERRQ(ierr);
  ierr = VecScale(final,a);CHKERRQ(ierr);

  VecAssemblyBegin(ini);
  VecAssemblyEnd(ini);

  VecAssemblyBegin(final);
  VecAssemblyEnd(final);

/*    VecView(final,PETSC_VIEWER_STDOUT_WORLD); */
  ierr = VecAXPY(final,-1.0,ini);CHKERRQ(ierr);
  ierr = VecNorm(final,NORM_1,&enorm);CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"  Error norm of |x - z|  = %e\n",enorm);CHKERRQ(ierr);
  fftw_destroy_plan(fplan);
  fftw_destroy_plan(bplan);
  fftw_free(in1); ierr = VecDestroy(&fin);CHKERRQ(ierr);
  fftw_free(out); ierr = VecDestroy(&fout);CHKERRQ(ierr);
  fftw_free(in2); ierr = VecDestroy(&fout1);CHKERRQ(ierr);

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


示例4: main

int main(int argc, char **argv)
{
  PetscInt    ierr;
  PetscSF     sf;
  Vec         A,Aout;
  Vec         B,Bout;
  PetscScalar *bufA;
  PetscScalar *bufAout;
  PetscScalar *bufB;
  PetscScalar *bufBout;
  PetscMPIInt rank, size;
  PetscInt    nroots, nleaves;
  PetscInt    i;
  PetscInt    *ilocal;
  PetscSFNode *iremote;

  ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
  ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);

  if (size != 2) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_USER, "Only coded for two MPI processes\n");

  ierr = PetscSFCreate(PETSC_COMM_WORLD,&sf);CHKERRQ(ierr);
  ierr = PetscSFSetFromOptions(sf);CHKERRQ(ierr);

  nleaves = 2;
  nroots = 1;
  ierr = PetscMalloc1(nleaves,&ilocal);CHKERRQ(ierr);

  for (i = 0; i<nleaves; i++) {
    ilocal[i] = i;
  }

  ierr = PetscMalloc1(nleaves,&iremote);CHKERRQ(ierr);
  if (rank == 0) {
    iremote[0].rank = 0;
    iremote[0].index = 0;
    iremote[1].rank = 1;
    iremote[1].index = 0;
  } else {
    iremote[0].rank = 1;
    iremote[0].index = 0;
    iremote[1].rank = 0;
    iremote[1].index = 0;
  }
  ierr = PetscSFSetGraph(sf,nroots,nleaves,ilocal,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);CHKERRQ(ierr);
  ierr = PetscSFSetUp(sf);CHKERRQ(ierr);
  ierr = PetscSFView(sf,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = VecSetSizes(A,2,PETSC_DETERMINE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(A);CHKERRQ(ierr);
  ierr = VecSetUp(A);CHKERRQ(ierr);

  ierr = VecDuplicate(A,&B);CHKERRQ(ierr);
  ierr = VecDuplicate(A,&Aout);CHKERRQ(ierr);
  ierr = VecDuplicate(A,&Bout);CHKERRQ(ierr);
  ierr = VecGetArray(A,&bufA);CHKERRQ(ierr);
  ierr = VecGetArray(B,&bufB);CHKERRQ(ierr);
  for (i=0; i<2; i++) {
    bufA[i] = (PetscScalar)rank;
    bufB[i] = (PetscScalar)(rank) + 10.0;
  }
  ierr = VecRestoreArray(A,&bufA);CHKERRQ(ierr);
  ierr = VecRestoreArray(B,&bufB);CHKERRQ(ierr);

  ierr = VecGetArrayRead(A,(const PetscScalar**)&bufA);CHKERRQ(ierr);
  ierr = VecGetArrayRead(B,(const PetscScalar**)&bufB);CHKERRQ(ierr);
  ierr = VecGetArray(Aout,&bufAout);CHKERRQ(ierr);
  ierr = VecGetArray(Bout,&bufBout);CHKERRQ(ierr);
  ierr = PetscSFBcastBegin(sf,MPIU_SCALAR,(const void*)bufA,(void *)bufAout);CHKERRQ(ierr);
  ierr = PetscSFBcastBegin(sf,MPIU_SCALAR,(const void*)bufB,(void *)bufBout);CHKERRQ(ierr);
  ierr = PetscSFBcastEnd(sf,MPIU_SCALAR,(const void*)bufA,(void *)bufAout);CHKERRQ(ierr);
  ierr = PetscSFBcastEnd(sf,MPIU_SCALAR,(const void*)bufB,(void *)bufBout);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(A,(const PetscScalar**)&bufA);CHKERRQ(ierr);
  ierr = VecRestoreArrayRead(B,(const PetscScalar**)&bufB);CHKERRQ(ierr);
  ierr = VecRestoreArray(Aout,&bufAout);CHKERRQ(ierr);
  ierr = VecRestoreArray(Bout,&bufBout);CHKERRQ(ierr);

  ierr = VecView(Aout,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecView(Bout,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
  ierr = VecDestroy(&A);CHKERRQ(ierr);
  ierr = VecDestroy(&B);CHKERRQ(ierr);
  ierr = VecDestroy(&Aout);CHKERRQ(ierr);
  ierr = VecDestroy(&Bout);CHKERRQ(ierr);
  ierr = PetscSFDestroy(&sf);CHKERRQ(ierr);

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


示例5: MatCholeskyFactorNumeric_SeqSBAIJ_4

/* Version for when blocks are 4 by 4  */
PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat C,Mat A,const MatFactorInfo *info)
{
  Mat_SeqSBAIJ   *a   = (Mat_SeqSBAIJ*)A->data,*b = (Mat_SeqSBAIJ*)C->data;
  IS             perm = b->row;
  PetscErrorCode ierr;
  const PetscInt *ai,*aj,*perm_ptr,mbs=a->mbs,*bi=b->i,*bj=b->j;
  PetscInt       i,j,*a2anew,k,k1,jmin,jmax,*jl,*il,vj,nexti,ili;
  MatScalar      *ba = b->a,*aa,*ap,*dk,*uik;
  MatScalar      *u,*diag,*rtmp,*rtmp_ptr;
  PetscBool      pivotinblocks = b->pivotinblocks;
  PetscReal      shift         = info->shiftamount;
  PetscBool      allowzeropivot,zeropivotdetected=PETSC_FALSE;

  PetscFunctionBegin;
  /* initialization */
  allowzeropivot = PetscNot(A->erroriffailure);
  ierr = PetscCalloc1(16*mbs,&rtmp);CHKERRQ(ierr);
  ierr = PetscMalloc2(mbs,&il,mbs,&jl);CHKERRQ(ierr);
  il[0] = 0;
  for (i=0; i<mbs; i++) jl[i] = mbs; 
  
  ierr = PetscMalloc2(16,&dk,16,&uik);CHKERRQ(ierr);
  ierr = ISGetIndices(perm,&perm_ptr);CHKERRQ(ierr);

  /* check permutation */
  if (!a->permute) {
    ai = a->i; aj = a->j; aa = a->a;
  } else {
    ai   = a->inew; aj = a->jnew;
    ierr = PetscMalloc1(16*ai[mbs],&aa);CHKERRQ(ierr);
    ierr = PetscMemcpy(aa,a->a,16*ai[mbs]*sizeof(MatScalar));CHKERRQ(ierr);
    ierr = PetscMalloc1(ai[mbs],&a2anew);CHKERRQ(ierr);
    ierr = PetscMemcpy(a2anew,a->a2anew,(ai[mbs])*sizeof(PetscInt));CHKERRQ(ierr);

    for (i=0; i<mbs; i++) {
      jmin = ai[i]; jmax = ai[i+1];
      for (j=jmin; j<jmax; j++) {
        while (a2anew[j] != j) {
          k = a2anew[j]; a2anew[j] = a2anew[k]; a2anew[k] = k;
          for (k1=0; k1<16; k1++) {
            dk[k1]      = aa[k*16+k1];
            aa[k*16+k1] = aa[j*16+k1];
            aa[j*16+k1] = dk[k1];
          }
        }
        /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
        if (i > aj[j]) {
          /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
          ap = aa + j*16;                     /* ptr to the beginning of j-th block of aa */
          for (k=0; k<16; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
          for (k=0; k<4; k++) {               /* j-th block of aa <- dk^T */
            for (k1=0; k1<4; k1++) *ap++ = dk[k + 4*k1];
          }
        }
      }
    }
    ierr = PetscFree(a2anew);CHKERRQ(ierr);
  }

  /* for each row k */
  for (k = 0; k<mbs; k++) {

    /*initialize k-th row with elements nonzero in row perm(k) of A */
    jmin = ai[perm_ptr[k]]; jmax = ai[perm_ptr[k]+1];
    if (jmin < jmax) {
      ap = aa + jmin*16;
      for (j = jmin; j < jmax; j++) {
        vj       = perm_ptr[aj[j]];   /* block col. index */
        rtmp_ptr = rtmp + vj*16;
        for (i=0; i<16; i++) *rtmp_ptr++ = *ap++;
      }
    }

    /* modify k-th row by adding in those rows i with U(i,k) != 0 */
    ierr = PetscMemcpy(dk,rtmp+k*16,16*sizeof(MatScalar));CHKERRQ(ierr);
    i    = jl[k]; /* first row to be added to k_th row  */

    while (i < mbs) {
      nexti = jl[i]; /* next row to be added to k_th row */

      /* compute multiplier */
      ili = il[i];  /* index of first nonzero element in U(i,k:bms-1) */

      /* uik = -inv(Di)*U_bar(i,k) */
      diag = ba + i*16;
      u    = ba + ili*16;

      uik[0] = -(diag[0]*u[0] + diag[4]*u[1] + diag[8]*u[2] + diag[12]*u[3]);
      uik[1] = -(diag[1]*u[0] + diag[5]*u[1] + diag[9]*u[2] + diag[13]*u[3]);
      uik[2] = -(diag[2]*u[0] + diag[6]*u[1] + diag[10]*u[2]+ diag[14]*u[3]);
      uik[3] = -(diag[3]*u[0] + diag[7]*u[1] + diag[11]*u[2]+ diag[15]*u[3]);

      uik[4] = -(diag[0]*u[4] + diag[4]*u[5] + diag[8]*u[6] + diag[12]*u[7]);
      uik[5] = -(diag[1]*u[4] + diag[5]*u[5] + diag[9]*u[6] + diag[13]*u[7]);
      uik[6] = -(diag[2]*u[4] + diag[6]*u[5] + diag[10]*u[6]+ diag[14]*u[7]);
      uik[7] = -(diag[3]*u[4] + diag[7]*u[5] + diag[11]*u[6]+ diag[15]*u[7]);

      uik[8] = -(diag[0]*u[8] + diag[4]*u[9] + diag[8]*u[10] + diag[12]*u[11]);
      uik[9] = -(diag[1]*u[8] + diag[5]*u[9] + diag[9]*u[10] + diag[13]*u[11]);
//.........这里部分代码省略.........
开发者ID:firedrakeproject,项目名称:petsc,代码行数:101,代码来源:sbaijfact6.c


示例6: main

int main(int argc,char **argv)
{
  TS                ts;         /* time integrator */
  TSAdapt           adapt;
  Vec               X;          /* solution vector */
  Mat               J;          /* Jacobian matrix */
  PetscInt          steps,maxsteps,ncells,xs,xm,i;
  PetscErrorCode    ierr;
  PetscReal         ftime,dt;
  char              chemfile[PETSC_MAX_PATH_LEN] = "chem.inp",thermofile[PETSC_MAX_PATH_LEN] = "therm.dat";
  struct _User      user;
  TSConvergedReason reason;
  PetscBool         showsolutions = PETSC_FALSE;
  char              **snames,*names;
  Vec               lambda;     /* used with TSAdjoint for sensitivities */

  ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Chemistry solver options","");CHKERRQ(ierr);
  ierr = PetscOptionsString("-chem","CHEMKIN input file","",chemfile,chemfile,sizeof(chemfile),NULL);CHKERRQ(ierr);
  ierr = PetscOptionsString("-thermo","NASA thermo input file","",thermofile,thermofile,sizeof(thermofile),NULL);CHKERRQ(ierr);
  user.pressure = 1.01325e5;    /* Pascal */
  ierr = PetscOptionsReal("-pressure","Pressure of reaction [Pa]","",user.pressure,&user.pressure,NULL);CHKERRQ(ierr);
  user.Tini   = 1550;
  ierr = PetscOptionsReal("-Tini","Initial temperature [K]","",user.Tini,&user.Tini,NULL);CHKERRQ(ierr);
  user.diffus = 100;
  ierr = PetscOptionsReal("-diffus","Diffusion constant","",user.diffus,&user.diffus,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsBool("-draw_solution","Plot the solution for each cell","",showsolutions,&showsolutions,NULL);CHKERRQ(ierr);
  user.diffusion = PETSC_TRUE;
  ierr = PetscOptionsBool("-diffusion","Have diffusion","",user.diffusion,&user.diffusion,NULL);CHKERRQ(ierr);
  user.reactions = PETSC_TRUE;
  ierr = PetscOptionsBool("-reactions","Have reactions","",user.reactions,&user.reactions,NULL);CHKERRQ(ierr);
  ierr = PetscOptionsEnd();CHKERRQ(ierr);

  ierr = TC_initChem(chemfile, thermofile, 0, 1.0);TCCHKERRQ(ierr);
  user.Nspec = TC_getNspec();
  user.Nreac = TC_getNreac();

  ierr    = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,-1,user.Nspec+1,1,NULL,&user.dm);CHKERRQ(ierr);
  ierr    = DMDAGetInfo(user.dm,NULL,&ncells,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
  user.dx = 1.0/ncells;  /* Set the coordinates of the cell centers; note final ghost cell is at x coordinate 1.0 */
  ierr    = DMDASetUniformCoordinates(user.dm,0.0,1.0,0.0,1.0,0.0,1.0);CHKERRQ(ierr);

  /* set the names of each field in the DMDA based on the species name */
  ierr = PetscMalloc1((user.Nspec+1)*LENGTHOFSPECNAME,&names);CHKERRQ(ierr);
  ierr = PetscStrcpy(names,"Temp");CHKERRQ(ierr);
  TC_getSnames(user.Nspec,names+LENGTHOFSPECNAME);CHKERRQ(ierr);
  ierr = PetscMalloc1((user.Nspec+2),&snames);CHKERRQ(ierr);
  for (i=0; i<user.Nspec+1; i++) snames[i] = names+i*LENGTHOFSPECNAME;
  snames[user.Nspec+1] = NULL;
  ierr = DMDASetFieldNames(user.dm,(const char * const *)snames);CHKERRQ(ierr);
  ierr = PetscFree(snames);CHKERRQ(ierr);
  ierr = PetscFree(names);CHKERRQ(ierr);


  ierr = DMCreateMatrix(user.dm,&J);CHKERRQ(ierr);
  ierr = DMCreateGlobalVector(user.dm,&X);CHKERRQ(ierr);

  ierr = PetscMalloc3(user.Nspec+1,&user.tchemwork,PetscSqr(user.Nspec+1),&user.Jdense,user.Nspec+1,&user.rows);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Create timestepping solver context
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
  ierr = TSSetDM(ts,user.dm);CHKERRQ(ierr);
  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);
  ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr);
  ierr = TSARKIMEXSetType(ts,TSARKIMEX4);CHKERRQ(ierr);
  ierr = TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);CHKERRQ(ierr);
  ierr = TSSetRHSJacobian(ts,J,J,FormRHSJacobian,&user);CHKERRQ(ierr);

  ftime    = 1.0;
  maxsteps = 10000;
  ierr     = TSSetDuration(ts,maxsteps,ftime);CHKERRQ(ierr);
  ierr     = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set initial conditions
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  ierr = FormInitialSolution(ts,X,&user);CHKERRQ(ierr);
  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
  dt   = 1e-10;                 /* Initial time step */
  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
  ierr = TSGetAdapt(ts,&adapt);CHKERRQ(ierr);
  ierr = TSAdaptSetStepLimits(adapt,1e-12,1e-4);CHKERRQ(ierr); /* Also available with -ts_adapt_dt_min/-ts_adapt_dt_max */
  ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr);            /* Retry step an unlimited number of times */


  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Pass information to graphical monitoring routine
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  if (showsolutions) {
    ierr = DMDAGetCorners(user.dm,&xs,NULL,NULL,&xm,NULL,NULL);CHKERRQ(ierr);
    for (i=xs;i<xs+xm;i++) {
      ierr = MonitorCell(ts,&user,i);CHKERRQ(ierr);
    }
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     Set runtime options
   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
//.........这里部分代码省略.........
开发者ID:tom-klotz,项目名称:petsc,代码行数:101,代码来源:extchemfield.c


示例7: PCGAMGCreateGraph

PetscErrorCode PCGAMGCreateGraph(const Mat Amat, Mat *a_Gmat)
{
  PetscErrorCode ierr;
  PetscInt       Istart,Iend,Ii,jj,kk,ncols,nloc,NN,MM,bs;
  PetscMPIInt    rank, size;
  MPI_Comm       comm;
  Mat            Gmat;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr);
  ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
  ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr);
  ierr = MatGetSize(Amat, &MM, &NN);CHKERRQ(ierr);
  ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr);
  nloc = (Iend-Istart)/bs;

#if defined PETSC_GAMG_USE_LOG
  ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
#endif
  if (bs > 1) {
    const PetscScalar *vals;
    const PetscInt    *idx;
    PetscInt          *d_nnz, *o_nnz;
    /* count nnz, there is sparcity in here so this might not be enough */
    ierr = PetscMalloc1(nloc, &d_nnz);CHKERRQ(ierr);
    ierr = PetscMalloc1(nloc, &o_nnz);CHKERRQ(ierr);
    for (Ii = Istart, jj = 0; Ii < Iend; Ii += bs, jj++) {
      d_nnz[jj] = 0;
      for (kk=0; kk<bs; kk++) {
        ierr = MatGetRow(Amat,Ii+kk,&ncols,0,0);CHKERRQ(ierr);
        if (ncols > d_nnz[jj]) {
          d_nnz[jj] = ncols; /* very pessimistic but could be too low in theory */
          o_nnz[jj] = ncols;
          if (d_nnz[jj] > nloc) d_nnz[jj] = nloc;
          if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc;
        }
        ierr = MatRestoreRow(Amat,Ii+kk,&ncols,0,0);CHKERRQ(ierr);
      }
    }

    /* get scalar copy (norms) of matrix -- AIJ specific!!! */
    ierr = MatCreateAIJ(comm, nloc, nloc, PETSC_DETERMINE, PETSC_DETERMINE,0, d_nnz, 0, o_nnz, &Gmat);CHKERRQ(ierr);

    ierr = PetscFree(d_nnz);CHKERRQ(ierr);
    ierr = PetscFree(o_nnz);CHKERRQ(ierr);

    for (Ii = Istart; Ii < Iend; Ii++) {
      PetscInt dest_row = Ii/bs;
      ierr = MatGetRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
      for (jj=0; jj<ncols; jj++) {
        PetscInt    dest_col = idx[jj]/bs;
        PetscScalar sv       = PetscAbs(PetscRealPart(vals[jj]));
        ierr = MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);CHKERRQ(ierr);
      }
      ierr = MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
  } else {
    /* just copy scalar matrix - abs() not taken here but scaled later */
    ierr = MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);CHKERRQ(ierr);
  }

#if defined PETSC_GAMG_USE_LOG
  ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr);
#endif

  *a_Gmat = Gmat;
  PetscFunctionReturn(0);
}
开发者ID:ZJLi2013,项目名称:petsc,代码行数:71,代码来源:tools.c


示例8: read

/*@C
   PetscBinaryRead - Reads from a binary file.

   Not Collective

   Input Parameters:
+  fd - the file descriptor
.  num  - the maximum number of items to read
-  type - the type of items to read (PETSC_INT, PETSC_REAL, PETSC_SCALAR, etc.)

   Output Parameters:
+  data - the buffer
-  count - the number of items read, optional



   Level: developer

   Notes:
   If count is not provided and the number of items read is less than
   the maximum number of items to read, then this routine errors.

   PetscBinaryRead() uses byte swapping to work on all machines; the files
   are written to file ALWAYS using big-endian ordering. On small-endian machines the numbers
   are converted to the small-endian format when they are read in from the file.
   When PETSc is ./configure with --with-64bit-indices the integers are written to the
   file as 64 bit integers, this means they can only be read back in when the option --with-64bit-indices
   is used.

   Concepts: files^reading binary
   Concepts: binary files^reading

.seealso: PetscBinaryWrite(), PetscBinaryOpen(), PetscBinaryClose(), PetscViewerBinaryGetDescriptor(), PetscBinarySynchronizedWrite(),
          PetscBinarySynchronizedRead(), PetscBinarySynchronizedSeek()
@*/
PetscErrorCode  PetscBinaryRead(int fd,void *data,PetscInt num,PetscInt *count,PetscDataType type)
{
  size_t            typesize, m = (size_t) num, n = 0, maxblock = 65536;
  char              *p = (char*)data;
#if defined(PETSC_USE_REAL___FLOAT128)
  PetscBool         readdouble = PETSC_FALSE;
  double            *pdouble;
#endif
  void              *ptmp = data;
  char              *fname = NULL;
  PetscErrorCode    ierr;

  PetscFunctionBegin;
  if (count) *count = 0;
  if (num < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to read a negative amount of data %D",num);
  if (!num) PetscFunctionReturn(0);

  if (type == PETSC_FUNCTION) {
    m     = 64;
    type  = PETSC_CHAR;
    fname = (char*)malloc(m*sizeof(char));
    p     = (char*)fname;
    ptmp  = (void*)fname;
    if (!fname) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MEM,"Cannot allocate space for function name");
  }
  if (type == PETSC_BIT_LOGICAL) m = PetscBTLength(m);

  ierr = PetscDataTypeGetSize(type,&typesize);CHKERRQ(ierr);

#if defined(PETSC_USE_REAL___FLOAT128)
  ierr = PetscOptionsGetBool(NULL,NULL,"-binary_read_double",&readdouble,NULL);CHKERRQ(ierr);
  /* If using __float128 precision we still read in doubles from file */
  if ((type == PETSC_REAL || type == PETSC_COMPLEX) && readdouble) {
    PetscInt cnt = num * ((type == PETSC_REAL) ? 1 : 2);
    ierr = PetscMalloc1(cnt,&pdouble);CHKERRQ(ierr);
    p = (char*)pdouble;
    typesize /= 2;
  }
#endif

  m *= typesize;

  while (m) {
    size_t len = (m < maxblock) ? m : maxblock;
    int    ret = (int)read(fd,p,len);
    if (ret < 0 && errno == EINTR) continue;
    if (!ret && len > 0) break; /* Proxy for EOF */
    if (ret < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Error reading from file, errno %d",errno);
    m -= ret;
    p += ret;
    n += ret;
  }
  if (m && !count) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Read past end of file");

  num = (PetscInt)(n/typesize); /* Should we require `n % typesize == 0` ? */
  if (count) *count = num;      /* TODO: This is most likely wrong for PETSC_BIT_LOGICAL */

#if defined(PETSC_USE_REAL___FLOAT128)
  if ((type == PETSC_REAL || type == PETSC_COMPLEX) && readdouble) {
    PetscInt  i, cnt = num * ((type == PETSC_REAL) ? 1 : 2);
    PetscReal *preal = (PetscReal*)data;
    if (!PetscBinaryBigEndian()) {ierr = PetscByteSwapDouble(pdouble,cnt);CHKERRQ(ierr);}
    for (i=0; i<cnt; i++) preal[i] = pdouble[i];
    ierr = PetscFree(pdouble);CHKERRQ(ierr);
    PetscFunctionReturn(0);
//.........这里部分代码省略.........
开发者ID:petsc,项目名称:petsc,代码行数:101,代码来源:sysio.c


示例9: PCBDDCNullSpaceAssembleCorrection


//.........这里部分代码省略.........
    ierr = VecScatterEnd(matis->rctx,nullvecs[k],pcis->vec1_N,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    if (basis_dofs) {
      ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
      ierr = VecScatterBegin(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
      ierr = VecScatterEnd(scatter_ctx,pcis->vec1_N,work1,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
      ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
      ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
      ierr = VecResetArray(work1);CHKERRQ(ierr);
      ierr = VecResetArray(work2);CHKERRQ(ierr);
    }
  }

  if (basis_dofs) {
    if (nnsp_has_cnst) {
      ierr = VecPlaceArray(work1,(const PetscScalar*)&basis_mat[k*basis_dofs]);CHKERRQ(ierr);
      ierr = VecSet(work1,one);CHKERRQ(ierr);
      ierr = VecPlaceArray(work2,(const PetscScalar*)&Kbasis_mat[k*basis_dofs]);CHKERRQ(ierr);
      ierr = MatMult(local_mat,work1,work2);CHKERRQ(ierr);
      ierr = VecResetArray(work1);CHKERRQ(ierr);
      ierr = VecResetArray(work2);CHKERRQ(ierr);
    }
    ierr = VecDestroy(&work1);CHKERRQ(ierr);
    ierr = VecDestroy(&work2);CHKERRQ(ierr);
    ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr);
    ierr = MatDenseRestoreArray(shell_ctx->basis_mat,&basis_mat);CHKERRQ(ierr);
    ierr = MatDenseRestoreArray(shell_ctx->Kbasis_mat,&Kbasis_mat);CHKERRQ(ierr);

    /* Assemble another Mat object in shell context */
    ierr = MatTransposeMatMult(shell_ctx->basis_mat,shell_ctx->Kbasis_mat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&small_mat);CHKERRQ(ierr);
    ierr = MatFactorInfoInitialize(&matinfo);CHKERRQ(ierr);
    ierr = ISCreateStride(PETSC_COMM_SELF,basis_size,0,1,&is_aux);CHKERRQ(ierr);
    ierr = MatLUFactor(small_mat,is_aux,is_aux,&matinfo);CHKERRQ(ierr);
    ierr = ISDestroy(&is_aux);CHKERRQ(ierr);
    ierr = PetscMalloc1(basis_size*basis_size,&array_mat);CHKERRQ(ierr);
    for (k=0;k<basis_size;k++) {
      ierr = VecSet(shell_ctx->work_small_1,zero);CHKERRQ(ierr);
      ierr = VecSetValue(shell_ctx->work_small_1,k,one,INSERT_VALUES);CHKERRQ(ierr);
      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);
开发者ID:plguhur,项目名称:petsc,代码行数:67,代码来源:bddcnullspace.c


示例10: PCBDDCNullSpaceAssembleCoarse

PetscErrorCode PCBDDCNullSpaceAssembleCoarse(PC pc, Mat coarse_mat, MatNullSpace* CoarseNullSpace)
{
  PC_BDDC        *pcbddc = (PC_BDDC*)pc->data;
  Mat_IS         *matis = (Mat_IS*)pc->pmat->data;
  MatNullSpace   tempCoarseNullSpace=NULL;
  const Vec      *nsp_vecs;
  Vec            *coarse_nsp_vecs,local_vec,local_primal_vec,wcoarse_vec,wcoarse_rhs;
  PetscInt       nsp_size,coarse_nsp_size,i;
  PetscBool      nsp_has_cnst;
  PetscReal      test_null;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  tempCoarseNullSpace = 0;
  coarse_nsp_size = 0;
  coarse_nsp_vecs = 0;
  ierr = MatNullSpaceGetVecs(pcbddc->NullSpace,&nsp_has_cnst,&nsp_size,&nsp_vecs);CHKERRQ(ierr);
  if (coarse_mat) {
    ierr = PetscMalloc1(nsp_size+1,&coarse_nsp_vecs);CHKERRQ(ierr);
    for (i=0;i<nsp_size+1;i++) {
      ierr = MatCreateVecs(coarse_mat,&coarse_nsp_vecs[i],NULL);CHKERRQ(ierr);
    }
    if (pcbddc->dbg_flag) {
      ierr = MatCreateVecs(coarse_mat,&wcoarse_vec,&wcoarse_rhs);CHKERRQ(ierr);
    }
  }
  ierr = MatCreateVecs(pcbddc->ConstraintMatrix,&local_vec,&local_primal_vec);CHKERRQ(ierr);
  if (nsp_has_cnst) {
    ierr = VecSet(local_vec,1.0);CHKERRQ(ierr);
    ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
    ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    if (coarse_mat) {
      PetscScalar *array_out;
      const PetscScalar *array_in;
      PetscInt lsize;
      if (pcbddc->dbg_flag) {
        PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
        ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr);
        ierr = VecNorm(wcoarse_rhs,NORM_INFINITY,&test_null);CHKERRQ(ierr);
        ierr = PetscViewerASCIIPrintf(dbg_viewer,"Constant coarse null space error % 1.14e\n",test_null);CHKERRQ(ierr);
        ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr);
      }
      ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr);
      ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
      ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
      ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr);
      ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
      ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
      coarse_nsp_size++;
    }
  }
  for (i=0;i<nsp_size;i++)  {
    ierr = VecScatterBegin(matis->rctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(matis->rctx,nsp_vecs[i],local_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = MatMult(pcbddc->ConstraintMatrix,local_vec,local_primal_vec);CHKERRQ(ierr);
    ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,local_primal_vec,pcbddc->coarse_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
    if (coarse_mat) {
      PetscScalar *array_out;
      const PetscScalar *array_in;
      PetscInt lsize;
      if (pcbddc->dbg_flag) {
        PetscViewer dbg_viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)coarse_mat));
        ierr = MatMult(coarse_mat,wcoarse_vec,wcoarse_rhs);CHKERRQ(ierr);
        ierr = VecNorm(wcoarse_rhs,NORM_2,&test_null);CHKERRQ(ierr);
        ierr = PetscViewerASCIIPrintf(dbg_viewer,"Vec %d coarse null space error % 1.14e\n",i,test_null);CHKERRQ(ierr);
        ierr = PetscViewerFlush(dbg_viewer);CHKERRQ(ierr);
      }
      ierr = VecGetLocalSize(pcbddc->coarse_vec,&lsize);CHKERRQ(ierr);
      ierr = VecGetArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
      ierr = VecGetArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
      ierr = PetscMemcpy(array_out,array_in,lsize*sizeof(PetscScalar));CHKERRQ(ierr);
      ierr = VecRestoreArray(coarse_nsp_vecs[coarse_nsp_size],&array_out);CHKERRQ(ierr);
      ierr = VecRestoreArrayRead(pcbddc->coarse_vec,&array_in);CHKERRQ(ierr);
      coarse_nsp_size++;
    }
  }
  if (coarse_nsp_size > 0) {
    ierr = PCBDDCOrthonormalizeVecs(coarse_nsp_size,coarse_nsp_vecs);CHKERRQ(ierr);
    ierr = MatNullSpaceCreate(PetscObjectComm((PetscObject)coarse_mat),PETSC_FALSE,coarse_nsp_size,coarse_nsp_vecs,&tempCoarseNullSpace);CHKERRQ(ierr);
    for (i=0;i<nsp_size+1;i++) {
      ierr = VecDestroy(&coarse_nsp_vecs[i]);CHKERRQ(ierr);
    }
  }
  if (coarse_mat) {
    ierr = PetscFree(coarse_nsp_vecs);CHKERRQ(ierr);
    if (pcbddc->dbg_flag) {
      ierr = VecDestroy(&wcoarse_vec);CHKERRQ(ierr);
      ierr = VecDestroy(&wcoarse_rhs);CHKERRQ(ierr);
    }
  }
  ierr = VecDestroy(&local_vec);CHKERRQ(ierr);
  ierr = VecDestroy(&local_primal_vec);CHKERRQ(ierr);
  *CoarseNullSpace = tempCoarseNullSpace;
  PetscFunctionReturn(0);
}
开发者ID:plguhur,项目名称:petsc,代码行数:97,代码来源:bddcnullspace.c


示例11: MatISGetMPIXAIJ_IS

PetscErrorCode MatISGetMPIXAIJ_IS(Mat mat, MatReuse reuse, Mat *M)
{
  Mat_IS         *matis = (Mat_IS*)(mat->data);
  Mat            local_mat;
  /* info on mat */
  PetscInt       bs,rows,cols,lrows,lcols;
  PetscInt       local_rows,local_cols;
  PetscBool      isdense,issbaij,isseqaij;
  PetscMPIInt    nsubdomains;
  /* values insertion */
  PetscScalar    *array;
  /* work */
  PetscErrorCode ierr;

  PetscFunctionBegin;
  /* get info from mat */
  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&nsubdomains);CHKERRQ(ierr);
  if (nsubdomains == 1) {
    if (reuse == MAT_INITIAL_MATRIX) {
      ierr = MatDuplicate(matis->A,MAT_COPY_VALUES,&(*M));CHKERRQ(ierr);
    } else {
      ierr = MatCopy(matis->A,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
    }
    PetscFunctionReturn(0);
  }
  ierr = MatGetSize(mat,&rows,&cols);CHKERRQ(ierr);
  ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr);
  ierr = MatGetLocalSize(mat,&lrows,&lcols);CHKERRQ(ierr);
  ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQAIJ,&isseqaij);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr);

  if (reuse == MAT_INITIAL_MATRIX) {
    MatType     new_mat_type;
    PetscBool   issbaij_red;

    /* determining new matrix type */
    ierr = MPIU_Allreduce(&issbaij,&issbaij_red,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
    if (issbaij_red) {
      new_mat_type = MATSBAIJ;
    } else {
      if (bs>1) {
        new_mat_type = MATBAIJ;
      } else {
        new_mat_type = MATAIJ;
      }
    }

    ierr = MatCreate(PetscObjectComm((PetscObject)mat),M);CHKERRQ(ierr);
    ierr = MatSetSizes(*M,lrows,lcols,rows,cols);CHKERRQ(ierr);
    ierr = MatSetBlockSize(*M,bs);CHKERRQ(ierr);
    ierr = MatSetType(*M,new_mat_type);CHKERRQ(ierr);
    ierr = MatISSetMPIXAIJPreallocation_Private(mat,*M,PETSC_FALSE);CHKERRQ(ierr);
  } else {
    PetscInt mbs,mrows,mcols,mlrows,mlcols;
    /* some checks */
    ierr = MatGetBlockSize(*M,&mbs);CHKERRQ(ierr);
    ierr = MatGetSize(*M,&mrows,&mcols);CHKERRQ(ierr);
    ierr = MatGetLocalSize(*M,&mlrows,&mlcols);CHKERRQ(ierr);
    if (mrows != rows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of rows (%d != %d)",rows,mrows);
    if (mcols != cols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of cols (%d != %d)",cols,mcols);
    if (mlrows != lrows) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local rows (%d != %d)",lrows,mlrows);
    if (mlcols != lcols) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong number of local cols (%d != %d)",lcols,mlcols);
    if (mbs != bs) SETERRQ2(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot reuse matrix. Wrong block size (%d != %d)",bs,mbs);
    ierr = MatZeroEntries(*M);CHKERRQ(ierr);
  }

  if (issbaij) {
    ierr = MatConvert(matis->A,MATSEQBAIJ,MAT_INITIAL_MATRIX,&local_mat);CHKERRQ(ierr);
  } else {
    ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr);
    local_mat = matis->A;
  }

  /* Set values */
  ierr = MatSetLocalToGlobalMapping(*M,mat->rmap->mapping,mat->cmap->mapping);CHKERRQ(ierr);
  if (isdense) { /* special case for dense local matrices */
    PetscInt i,*dummy_rows,*dummy_cols;

    if (local_rows != local_cols) {
      ierr = PetscMalloc2(local_rows,&dummy_rows,local_cols,&dummy_cols);CHKERRQ(ierr);
      for (i=0;i<local_rows;i++) dummy_rows[i] = i;
      for (i=0;i<local_cols;i++) dummy_cols[i] = i;
    } else {
      ierr = PetscMalloc1(local_rows,&dummy_rows);CHKERRQ(ierr);
      for (i=0;i<local_rows;i++) dummy_rows[i] = i;
      dummy_cols = dummy_rows;
    }
    ierr = MatSetOption(*M,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr);
    ierr = MatDenseGetArray(local_mat,&array);CHKERRQ(ierr);
    ierr = MatSetValuesLocal(*M,local_rows,dummy_rows,local_cols,dummy_cols,array,ADD_VALUES);CHKERRQ(ierr);
    ierr = MatDenseRestoreArray(local_mat,&array);CHKERRQ(ierr);
    if (dummy_rows != dummy_cols) {
      ierr = PetscFree2(dummy_rows,dummy_cols);CHKERRQ(ierr);
    } else {
      ierr = PetscFree(dummy_rows);CHKERRQ(ierr);
    }
  } else if (isseqaij) {
    PetscInt  i,nvtxs,*xadj,*adjncy;
//.........这里部分代码省略.........
开发者ID:masa-ito,项目名称:PETScToPoisson,代码行数:101,代码来源:matis.c


示例12: MatISSetMPIXAIJPreallocation_Private

PETSC_EXTERN PetscErrorCode MatISSetMPIXAIJPreallocation_Private(Mat A, Mat B, PetscBool maxreduce)
{
  Mat_IS          *matis = (Mat_IS*)(A->data);
  PetscInt        *my_dnz,*my_onz,*dnz,*onz,*mat_ranges,*row_ownership;
  const PetscInt  *global_indices_r,*global_indices_c;
  PetscInt        i,j,bs,rows,cols;
  PetscInt        lrows,lcols;
  PetscInt        local_rows,local_cols;
  PetscMPIInt     nsubdomains;
  PetscBool       isdense,issbaij;
  PetscErrorCode  ierr;

  PetscFunctionBegin;
  ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&nsubdomains);CHKERRQ(ierr);
  ierr = MatGetSize(A,&rows,&cols);CHKERRQ(ierr);
  ierr = MatGetBlockSize(A,&bs);CHKERRQ(ierr);
  ierr = MatGetSize(matis->A,&local_rows,&local_cols);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQDENSE,&isdense);CHKERRQ(ierr);
  ierr = PetscObjectTypeCompare((PetscObject)matis->A,MATSEQSBAIJ,&issbaij);CHKERRQ(ierr);
  ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_r);CHKERRQ(ierr);
  if (A->rmap->mapping != A->cmap->mapping) {
    ierr = ISLocalToGlobalMappingGetIndices(A->rmap->mapping,&global_indices_c);CHKERRQ(ierr);
  } else {
    global_indices_c = global_indices_r;
  }

  if (issbaij) {
    ierr = MatGetRowUpperTriangular(matis->A);CHKERRQ(ierr);
  }
  /*
     An SF reduce is needed to sum up properly on shared rows.
     Note that generally preallocation is not exact, since it overestimates nonzeros
  */
  if (!matis->sf) { /* setup SF if not yet created and allocate rootdata and leafdata */
    ierr = MatISComputeSF_Private(A);CHKERRQ(ierr);
  }
  ierr = MatGetLocalSize(A,&lrows,&lcols);CHKERRQ(ierr);
  ierr = MatPreallocateInitialize(PetscObjectComm((PetscObject)A),lrows,lcols,dnz,onz);CHKERRQ(ierr);
  /* All processes need to compute entire row ownership */
  ierr = PetscMalloc1(rows,&row_ownership);CHKERRQ(ierr);
  ierr = MatGe 

鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
C++ PetscMalloc2函数代码示例发布时间:2022-05-30
下一篇:
C++ PetscMalloc函数代码示例发布时间:2022-05-30
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap