Commits

Hong Zhang  committed 9bfd627

bugfix of cholesky factorization with reordering for seqaij matrix

Hg-commit: b4b752b5a31797ea5a849cd6dbe758f544a8d7b7

  • Participants
  • Parent commits 70097a7

Comments (0)

Files changed (3)

File src/ksp/pc/impls/factor/cholesky/cholesky.c

     MatInfo info;
     if (!pc->setupcalled) {
       ierr = MatGetOrdering(pc->pmat,dir->ordering,&dir->row,&dir->col);CHKERRQ(ierr);
-      if (dir->col && (dir->row != dir->col)) {  /* only use row ordering for SBAIJ */
-        ierr = ISDestroy(dir->col);CHKERRQ(ierr); 
-        dir->col=0; 
-      }
+      /* check if dir->row == dir->col */
+      ierr = ISEqual(dir->row,dir->col,&flg);CHKERRQ(ierr);
+      if (!flg) SETERRQ(PETSC_ERR_ARG_INCOMP,"row and column permutations must equal");
+      ierr = ISDestroy(dir->col);CHKERRQ(ierr); /* only pass one ordering into CholeskyFactor */
+      dir->col=0; 
+
       ierr = PetscOptionsHasName(pc->prefix,"-pc_factor_nonzeros_along_diagonal",&flg);CHKERRQ(ierr);
       if (flg) {
         PetscReal tol = 1.e-10;

File src/mat/impls/aij/seq/aijfact.c

   Mat            C = *B;
   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
   Mat_SeqSBAIJ   *b=(Mat_SeqSBAIJ*)C->data;
-  IS             ip=b->row;
+  IS             ip=b->row,iip = b->icol;
   PetscErrorCode ierr;
-  PetscInt       *rip,i,j,mbs=A->rmap.n,*bi=b->i,*bj=b->j,*bcol;
+  PetscInt       *rip,*riip,i,j,mbs=A->rmap.n,*bi=b->i,*bj=b->j,*bcol;
   PetscInt       *ai=a->i,*aj=a->j;
   PetscInt       k,jmin,jmax,*jl,*il,col,nexti,ili,nz;
   MatScalar      *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi;
   zeropivot = info->zeropivot; 
 
   ierr  = ISGetIndices(ip,&rip);CHKERRQ(ierr);
+  ierr  = ISGetIndices(iip,&riip);CHKERRQ(ierr);
   
   /* initialization */
   nz   = (2*mbs+1)*sizeof(PetscInt)+mbs*sizeof(MatScalar);
       /* initialize k-th row by the perm[k]-th row of A */
       jmin = ai[rip[k]]; jmax = ai[rip[k]+1];
       for (j = jmin; j < jmax; j++){
-        col = rip[aj[j]];
+        col = riip[aj[j]];
         if (col >= k){ /* only take upper triangular entry */
           rtmp[col] = aa[j];
           *bval++  = 0.0; /* for in-place factorization */
   ierr = PetscFree(il);CHKERRQ(ierr);
 
   ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(iip,&riip);CHKERRQ(ierr);
   C->factor       = FACTOR_CHOLESKY; 
   C->assembled    = PETSC_TRUE; 
   C->preallocated = PETSC_TRUE;
 
   PetscFunctionBegin;
   /* check whether perm is the identity mapping */
-  ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);
+  ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr);  
+  ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr);
+  ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr);  
   ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr);
 
-  if (!perm_identity){
-    /* check if perm is symmetric! */
-    ierr = ISInvertPermutation(perm,PETSC_DECIDE,&iperm);CHKERRQ(ierr);  
-    ierr = ISGetIndices(iperm,&riip);CHKERRQ(ierr);
-    for (i=0; i<am; i++) {
-      if (rip[i] != riip[i]) SETERRQ(PETSC_ERR_ARG_INCOMP,"Non-symmetric permutation, must use symmetric permutation");
-    }
-    ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr);
-    ierr = ISDestroy(iperm);CHKERRQ(ierr);
-  } 
-
   /* initialization */
   ierr  = PetscMalloc((am+1)*sizeof(PetscInt),&ui);CHKERRQ(ierr);
   ui[0] = 0; 
     /* initialize lnk by the column indices of row rip[k] of A */
     nzk   = 0;
     ncols = ai[rip[k]+1] - ai[rip[k]]; 
+    if (!ncols) SETERRQ(PETSC_ERR_MAT_CH_ZRPVT,"Empty row in matrix");
     ncols_upper = 0;
     for (j=0; j<ncols; j++){
-      i = rip[*(aj + ai[rip[k]] + j)];
+      i = riip[*(aj + ai[rip[k]] + j)];  
       if (i >= k){ /* only take upper triangular entry */
         cols[ncols_upper] = i;
         ncols_upper++;
 #endif
 
   ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr);
+  ierr = ISRestoreIndices(iperm,&riip);CHKERRQ(ierr);
   ierr = PetscFree(jl);CHKERRQ(ierr);
 
   /* destroy list of free space and other temporary array(s) */
   b->ilen = 0;
   b->imax = 0;
   b->row  = perm;
-  b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
+  b->col  = perm;
   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 
-  b->icol = perm;
   ierr    = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 
+  b->icol = iperm;
+  b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */
   ierr    = PetscMalloc((am+1)*sizeof(PetscScalar),&b->solve_work);CHKERRQ(ierr);
   ierr    = PetscLogObjectMemory(B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr);
   b->maxnz = b->nz = ui[am];

File src/mat/impls/sbaij/seq/sbaij.c

   PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap.N,a->nz);
 #endif
   ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr);
-  if (a->row) {
-    ierr = ISDestroy(a->row);CHKERRQ(ierr);
-  }
+  if (a->row) {ierr = ISDestroy(a->row);CHKERRQ(ierr);}
+  if (a->col){ierr = ISDestroy(a->col);CHKERRQ(ierr);}
+  if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);}
   ierr = PetscFree(a->diag);CHKERRQ(ierr);
   ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr);
-  if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);}
   ierr = PetscFree(a->solve_work);CHKERRQ(ierr);
   ierr = PetscFree(a->relax_work);CHKERRQ(ierr);
   ierr = PetscFree(a->solves_work);CHKERRQ(ierr);