<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">Stefano :<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div dir="ltr">There's a bug in MatGetSubMatrix; the case MAT_REUSE_MATRIX is not covered by the interface when  isrow == iscol and isrow is a stride 1 IS representing the whole set of rows. A tentative patch is attached. </div></blockquote><div><br></div><div>Actually, current code is correct. The  case MAT_REUSE_MATRIX is covered:</div><div><div>          if (n == rend-rstart) {</div><div>            /* special case grabbing all rows; NEED to do a global reduction to make sure all processes are doing this */</div><div>            if (cll == MAT_INITIAL_MATRIX) {</div><div>              *newmat = mat;</div><div>              ierr    = PetscObjectReference((PetscObject)mat);CHKERRQ(ierr);</div><div>            }</div><div>            PetscFunctionReturn(0);</div><div>          }</div></div><div><br></div><div>Since 'newmat' is same as 'mat' in this case, we do NOT create a newmat, instead, assign it to the address of 'mat'  in the case MAT_INITIAL_MATRIX, i.e., </div><div>'*newmat' = mat are the same matrix. </div><div><br></div><div>Whenever 'mat' is changed, so does '*newmat' , MatGetSubMatrix(C,isrow,NULL,MAT_REUSE_MATRIX,&A) does not need to do anything.</div><div><br></div><div>I created a test ex181.c, attached below. </div><div><br></div><div>Hong</div></div></div></div>