In some cases after certain rows/columns have been removed due to active constraints certain OTHER rows end up being identically zero. This is some hacky code to remove those rows (since they screw up the linear solve). Just leave the code.
Barry On Dec 21, 2011, at 10:24 PM, Jed Brown wrote: > On Wed, Dec 21, 2011 at 19:16, Barry Smith <bsmith at mcs.anl.gov> wrote: > You could do this. Sounds reasonable. > > Done: > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/918daef122af > > > Now what is going on with this (apparent) debugging code virs.c:424 and > virsaug.c:1183 (http://petsc.cs.iit.edu/petsc/petsc-dev/rev/31c8ce): > > ierr = MatFindNonzeroRows(jac_inact_inact,&keptrows);CHKERRQ(ierr); > if (0 && keptrows) { > PetscInt cnt,*nrows,k; > const PetscInt *krows,*inact; > PetscInt rstart=jac_inact_inact->rmap->rstart; > > ierr = MatDestroy(&jac_inact_inact);CHKERRQ(ierr); > ierr = ISDestroy(&IS_act);CHKERRQ(ierr); > > ierr = ISGetLocalSize(keptrows,&cnt);CHKERRQ(ierr); > ierr = ISGetIndices(keptrows,&krows);CHKERRQ(ierr); > ierr = ISGetIndices(IS_inact,&inact);CHKERRQ(ierr); > ierr = PetscMalloc(cnt*sizeof(PetscInt),&nrows);CHKERRQ(ierr); > for (k=0; k<cnt; k++) { > nrows[k] = inact[krows[k]-rstart]; > } > ierr = ISRestoreIndices(keptrows,&krows);CHKERRQ(ierr); > ierr = ISRestoreIndices(IS_inact,&inact);CHKERRQ(ierr); > ierr = ISDestroy(&keptrows);CHKERRQ(ierr); > ierr = ISDestroy(&IS_inact);CHKERRQ(ierr); > > ierr = > ISCreateGeneral(PETSC_COMM_WORLD,cnt,nrows,PETSC_OWN_POINTER,&IS_inact);CHKERRQ(ierr); > ierr = > ISComplement(IS_inact,F->map->rstart,F->map->rend,&IS_act);CHKERRQ(ierr); > ierr = > MatGetSubMatrix(snes->jacobian,IS_inact,IS_inact,MAT_INITIAL_MATRIX,&jac_inact_inact);CHKERRQ(ierr); > } > ierr = DMSetVI(snes->dm,IS_inact);CHKERRQ(ierr); > /* remove later */ > > > This (a) leaks keptrows and (b) doesn't do anything useful because of the if > (0 && ...). What's the idea here?
