Showing results for

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted
##

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.

Sören_Plönnigs

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-08-2010
03:14 PM

46 Views

Some questions on pardiso

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

6 Replies

Highlighted
##

Hi Soeren,

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-08-2010
07:56 PM

46 Views

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

Highlighted
##

Hi Alexander,

Sören_Plönnigs

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2010
10:06 AM

46 Views

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

Highlighted
##

Hi Soeren,

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2010
10:56 AM

46 Views

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

Highlighted
##

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

Sören_Plönnigs

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2010
03:05 PM

46 Views

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

Highlighted
##

Hi Soeren,

Alexander_K_Intel2

Employee

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-11-2010
09:18 PM

46 Views

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

Highlighted
##

Sören_Plönnigs

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

08-12-2010
03:12 AM

46 Views

Hey Alex,

you are the man! With calloc everything works now.

Regards,

Soeren

For more complete information about compiler optimizations, see our Optimization Notice.