[petsc-users] Incomplete Cholesky Factorization in petsc

[Wenlong Gong]

Thanks, I'll look into it. -Wendy

On Tue, Dec 5, 2017 at 12:10 PM, Hong <hzhang at mcs.anl.gov> wrote:

> Wendy,
>
>>
>> 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?
>>
>
> You can look at MatSolve_SeqSBAIJ_1_NaturalOrdering() in
> petsc/src/mat/impls/sbaij/seq/sbaijfact2.c
>
>   /* solve U^T*D*y = b by forward substitution */
>   ierr = PetscMemcpy(x,b,mbs*sizeof(PetscScalar));CHKERRQ(ierr);
>   for (i=0; i<mbs; i++) {
>     v  = aa + ai[i];
>     vj = aj + ai[i];
>     xi = x[i];
>     nz = ai[i+1] - ai[i] - 1; /* exclude diag[i] */
>     for (j=0; j<nz; j++) x[vj[j]] += v[j]* xi;
>     x[i] = xi*v[nz];  /* v[nz] = aa[diag[i]] = 1/D(i) */
>   }
>
> i -- row index
> vj[j] -- col index
> v[j] -- matrix value
>
> Note: the input of this routine is a matrix icc/cholesky factor.
>
> Hong
>
>>
>> -Wendy
>>
>> On Tue, Dec 5, 2017 at 9:50 AM, Hong <hzhang at mcs.anl.gov> wrote:
>>
>>> Wendy,
>>> Petsc factor matrix as
>>> A = U^T * U, U: upper triangular matrix.
>>> When A is stored in seqaij format, we store U as
>>> /*
>>>    icc() under revised new data structure.
>>>    Factored arrays bj and ba are stored as
>>>      U(0,:),...,U(i,:),U(n-1,:)
>>>
>>>    ui=fact->i is an array of size n+1, in which
>>>    ui+
>>>      ui[i]:  points to 1st entry of U(i,:),i=0,...,n-1
>>>      ui[n]:  points to U(n-1,n-1)+1
>>>
>>>   udiag=fact->diag is an array of size n,in which
>>>      udiag[i]: points to diagonal of U(i,:), i=0,...,n-1
>>>
>>>    U(i,:) contains udiag[i] as its last entry, i.e.,
>>>     U(i,:) = (u[i,i+1],...,u[i,n-1],diag[i])
>>> */
>>>
>>> see petsc/src/mat/impls/aij/seq/aijfact.c: MatICCFactorSymbolic_SeqAIJ()
>>>
>>> Petsc matrix factors only work for Petsc MatSolve().
>>>
>>> Hong
>>>
>>> On Tue, Dec 5, 2017 at 9:32 AM, Wenlong Gong <wenlonggong at gmail.com>
>>> wrote:
>>>
>>>> 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.
>>>> 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().
>>>>
>>>> -Wendy
>>>>
>>>>
>>>> On Mon, Dec 4, 2017 at 10:07 PM, Smith, Barry F. <bsmith at mcs.anl.gov>
>>>> wrote:
>>>>
>>>>>
>>>>>   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.
>>>>>
>>>>>   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.
>>>>>
>>>>>     Barry
>>>>>
>>>>>
>>>>> > On Dec 4, 2017, at 9:07 PM, Wenlong Gong <wenlonggong at gmail.com>
>>>>> wrote:
>>>>> >
>>>>> > Hello,
>>>>> >
>>>>> > 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.
>>>>> >
>>>>> > 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!
>>>>> >
>>>>> > Best regards,
>>>>> > Wendy
>>>>> >
>>>>> > <repex.c>
>>>>>
>>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20171205/7cfca6e9/attachment.html>