<div dir="ltr">Wendy,<br><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div>Is there a function to get row/col indices and non-zeros vector of the factored matrix , instead of extracting the elements one by one?</div></div></blockquote><div> </div><div>You can look at MatSolve_SeqSBAIJ_1_NaturalOrdering() in</div><div>petsc/src/mat/impls/sbaij/seq/sbaijfact2.c </div><div><br></div><div><div> /* solve U^T*D*y = b by forward substitution */</div><div> ierr = PetscMemcpy(x,b,mbs*sizeof(PetscScalar));CHKERRQ(ierr);</div><div> for (i=0; i<mbs; i++) {</div><div> v = aa + ai[i];</div><div> vj = aj + ai[i];</div><div> xi = x[i];</div><div> nz = ai[i+1] - ai[i] - 1; /* exclude diag[i] */</div><div> for (j=0; j<nz; j++) x[vj[j]] += v[j]* xi;</div><div> x[i] = xi*v[nz]; /* v[nz] = aa[diag[i]] = 1/D(i) */</div><div> }</div></div><div><br></div><div>i -- row index</div><div>vj[j] -- col index</div><div>v[j] -- matrix value</div><div><br></div><div>Note: the input of this routine is a matrix icc/cholesky factor.</div><div><br></div><div>Hong</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div><br></div><div>-Wendy </div></div><div class="gmail-HOEnZb"><div class="gmail-h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Dec 5, 2017 at 9:50 AM, Hong <span dir="ltr"><<a href="mailto:hzhang@mcs.anl.gov" target="_blank">hzhang@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Wendy,<div>Petsc factor matrix as</div><div>A = U^T * U, U: upper triangular matrix.</div><div>When A is stored in seqaij format, we store U as</div><div><div>/*</div><div> icc() under revised new data structure.</div><div> Factored arrays bj and ba are stored as</div><div> U(0,:),...,U(i,:),U(n-1,:)</div><div><br></div><div> ui=fact->i is an array of size n+1, in which</div><div> ui+</div><div> ui[i]: points to 1st entry of U(i,:),i=0,...,n-1</div><div> ui[n]: points to U(n-1,n-1)+1</div><div><br></div><div> udiag=fact->diag is an array of size n,in which</div><div> udiag[i]: points to diagonal of U(i,:), i=0,...,n-1</div><div><br></div><div> U(i,:) contains udiag[i] as its last entry, i.e.,</div><div> U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i]<wbr>)</div><div>*/</div><div><br></div><div>see petsc/src/mat/impls/aij/se<wbr>q/aijfact.c: MatICCFactorSymbo<wbr>lic_SeqAIJ()</div><div><br></div><div>Petsc matrix factors only work for Petsc MatSolve(). </div><span class="gmail-m_-3942771939760444223HOEnZb"><font color="#888888"><div><br></div><div>Hong</div></font></span><div><div class="gmail-m_-3942771939760444223h5"><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Dec 5, 2017 at 9:32 AM, Wenlong Gong <span dir="ltr"><<a href="mailto:wenlonggong@gmail.com" target="_blank">wenlonggong@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Barry, Thank you for the reply! It makes sense If two IC() function returns slightly different factored matrix. But the IC factor matrix from PETSc is not even close to what I get from ichol in Matlab. <div>Ultimately I need this incomplete cholesky factorization function wrapped in R once it gives correct result. Not sure if I used the easiest way to get IC so I also seek for help on getting a concise version of IC().</div><div><br></div><div>-Wendy</div><div><br></div></div><div class="gmail-m_-3942771939760444223m_3016676108138644163gmail-HOEnZb"><div class="gmail-m_-3942771939760444223m_3016676108138644163gmail-h5"><div class="gmail_extra"><br><div class="gmail_quote">On Mon, Dec 4, 2017 at 10:07 PM, Smith, Barry F. <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
I'm not sure what your goal is. In general two different IC(0) codes might produce slightly different factorizations based on implementation details even if they claim to run the same algorithm so I don't think there is a reason to try to compare the factors they produce.<br>
<br>
In addition PETSc, as well as most IC/ILU/LU codes store the factored matrices in a "non-conventional" form that is optimized for performance so it is not easy to just pull out the "factors" to look at them or compare them.<br>
<br>
Barry<br>
<div><div class="gmail-m_-3942771939760444223m_3016676108138644163gmail-m_1166779876073856270h5"><br>
<br>
> On Dec 4, 2017, at 9:07 PM, Wenlong Gong <<a href="mailto:wenlonggong@gmail.com" target="_blank">wenlonggong@gmail.com</a>> wrote:<br>
><br>
> Hello,<br>
><br>
> I'm trying to use the Incomplete Cholesky Factorization for a sparse matrix in petsc. I started with a 10*10 matrix and used ksp and pc in order to get the ICC(0) factor, with no fill-in, natural ordering. However, the returned factor matrix does not match with the answer I got from matlab ichol() function.<br>
><br>
> The code with the hard-coded data is attached here. I would appreciate if anyone could help check if I did anything wrong.Please let me know if there is easier way to get this incomplete cholesky factor. Thanks!<br>
><br>
> Best regards,<br>
> Wendy<br>
><br>
</div></div>> <repex.c><br>
<br>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div></div></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br></div></div>