Some questions on pardiso

[Sören_Plönnigs]

I use the pardiso lib to solve the equation systems in my fem program. In my program I try to compare the direct algorithm with CG and CGS. My matrices are spd.

My question is if it is possible to use the CG/CGS solver without preconditioning? According to the manual, I tried to solve my system of equations by first running phase 11 followed by phase 33. In my eyes this should avoid the factorization. Why does the solver report the same number of nonzeros for matrix L as if I use phase 11 and 23 for solution? Shouldn't there be no matrix L in the first case? The solution takes two iterations for CG and only one iteration for CGS in both cases, which is an indication to me, that preconditioning is used in both cases.

Another Question I have is if the direct solver (phase 11 and 33) uses the same factorization for L and U then the CG/CGS solver? He reports only nonzeros for L, U has only one nonzero. Is this normal?

Regards,

Soeren

[Alexander_K_Intel2]

Hi Soeren,

Firstly, what mtype you use in you test? if you use mtype = 2 (spd matrices) and use iterative refinements that CG with preconditioner used as iterative method in other case - CGS with preconditioner. PARDISO use LU as preconditioner so if you want to use CG without any precontioner you could use cg subroutines from MKL (Intel Math Kernel Library ReferenceManual, Chapter 8, Itterative Sparce Solvers). About second qeustion: If you work with mtype = 2 then your LU decomposition provide U =L^T so number of nonzero elements in U as equal number of nonxero in L.

With best regards,

Alexander Kalinkin

[Sören_Plönnigs]

Hi Alexander,

I use mtype=2.

I now implemented a routine without preconditioning with iterative sparse solvers. but there is one problem, when the matrix gets bigger than 500000x500000 I get a segfault from the solver. Is there a maximum size for the solver? The system is 64bit and has 16GB of RAM, so this shouldn't be the problem.

And back to the direct Pardiso Solver, I also thought that the LU factorization should be LLT and so number(L)=number(U). In fact I get the following output from Pardiso:

number of nonzeros in L 18479

number of nonzeros in U 1

This seems a little bit strange to me, as the result I get from the solver is correct.

I also have one more question to the iterative Pardiso solvers, is the factorization that takes place complete, or is the fill-in limited? Is there a way to influence how much fill-in is allowed?

Regards,

Soeren

[Alexander_K_Intel2]

Hi Soeren,

Situation with iss is really strange, could you provide compiling line and testcase to check it?

What about PARDISO: number of nonzero in L/U is statistical information and print both LL^T and LU decomposition. So if you have LL^T decomposition the number of nonzero elements in U isunnecessary information. The PARDISO provide full LU decomposition of initial matrix after reorder, the incomplete factorization realized in MKL too, the functions described inIntel Math Kernel Library ReferenceManual, Chapter 8, Itterative Sparce Solvers,Preconditioners based on Incomplete LU Factorization Technique.

With best regards,

Alexander Kalinkin

[Sören_Plönnigs]

At the moment I only can give you the command line, which is

icpc -DMKL_LP64 -L /libmkl_solver_lp64.a -Wl,--start-group -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -Wl,--end-group -openmp -parallel -xSSSE3 -g -O0 -lpthread

providing an test case is not that easy, because i have no textfile of the matrix, it gets generated at an earlier time in the program. maybe you have already an idea when you see my commandline. by the way, I tried to run the program in gdb, which told me that the segfault comes from dcg_init:

0x000000000040e32c in iss_cg (A=@0x7fffdd3850f0, x=@0x7fffdd385180, b=@0x7fffdd385190) at solver/iterative.cpp:26

26 dcg_init(&n,x.values(),b.values(),&rci_request,ipar,dpar,tmp);

x and b are vector classes, the method value() gives access to the value array.

here's the part of my code that initializes the solver:

int ipar[128];

double dpar[128];

for (int i=0;i<128;i++) {

ipar = 0;

dpar = 0.0;

}

double tmp[4*n];

for (int i=0;i<4*n;i++)

tmp = 0;

int rci_request;

dcg_init(&n,x.values(),b.values(),&rci_request,ipar,dpar,tmp);

regards,

Soeren

[Alexander_K_Intel2]

Hi Soeren,

It's seem you have problem with lenght of statical array tmp. If you want to work with huge matrices try allocate array tmp, a, solution and so on.

With best regards,

Alexander Kalinkin

[Sören_Plönnigs]

Hey Alex,

you are the man! With calloc everything works now.

Regards,

Soeren