[petsc-users] Incomplete Cholesky Factorization in petsc

[Hong]

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/b9cab2cf/attachment.html>