topic Hi, in Intel® oneAPI Math Kernel Library

Conjugate Gradient source Code Example

Juan_Diego_E_ — Fri, 06 Feb 2015 23:50:20 GMT

Good afternoon,

I would like to have some help in the process of linking MKL with Visual Studio in Fortran Language. I am trying to run the example "cg_jacobi_precon" that comes in the Intel Folder, but I have not succeed. I am actually using F90. Can I have any advice, please, on which source codes do I need to include in my project as well as header files, and if I have to set any address in the Fortran and/or Linker properties of the project?

Please find below the actual example (I guess it is in F77) from Intel.

Thank you.

Regards,

JD.

!*******************************************************************************
!                              INTEL CONFIDENTIAL
!   Copyright(C) 2005-2010 Intel Corporation. All Rights Reserved.
!   The source code contained  or  described herein and all documents related to
!   the source code ("Material") are owned by Intel Corporation or its suppliers
!   or licensors.  Title to the  Material remains with  Intel Corporation or its
!   suppliers and licensors. The Material contains trade secrets and proprietary
!   and  confidential  information of  Intel or its suppliers and licensors. The
!   Material  is  protected  by  worldwide  copyright  and trade secret laws and
!   treaty  provisions. No part of the Material may be used, copied, reproduced,
!   modified, published, uploaded, posted, transmitted, distributed or disclosed
!   in any way without Intel's prior express written permission.
!   No license  under any  patent, copyright, trade secret or other intellectual
!   property right is granted to or conferred upon you by disclosure or delivery
!   of the Materials,  either expressly, by implication, inducement, estoppel or
!   otherwise.  Any  license  under  such  intellectual property  rights must be
!   express and approved by Intel in writing.
!
!*******************************************************************************
!  Content: Intel MKL RCI (P)CG Fortran-77 example
!
!*******************************************************************************

!---------------------------------------------------------------------------
!  Example program for solving symmetric positive definite system of equations.
!  Full case: full functionality of RCI (P)CG is used.
!---------------------------------------------------------------------------
      PROGRAM rci_pcg_f77_test_3
      IMPLICIT NONE
      INCLUDE 'mkl_rci.fi'
!---------------------------------------------------------------------------
! Define arrays for the upper triangle of the coefficient matrix and
! preconditioner as well as an array for rhs vector
! Compressed sparse row storage is used for sparse representation
!---------------------------------------------------------------------------
      INTEGER N, RCI_request, itercount, expected_itercount, i
      PARAMETER (N=8)
      PARAMETER (expected_itercount=8)
      DOUBLE PRECISION  rhs(N)
      INTEGER IA(9)
      INTEGER JA(18)
      DOUBLE PRECISION A(18)
! Fill all arrays containing matrix data.
      DATA IA /1,5,8,10,12,15,17,18,19/
      DATA JA
     1 /1,  3,     6,7,
     2    2,3,   5,
     3      3,         8,
     4         4,    7,
     5           5,6,7,
     6             6,  8,
     7               7,
     8                 8/
      DATA A
     1 /7.D0,      1.D0,           2.D0, 7.D0,
     2       -4.D0, 8.D0,     2.D0,
     3             1.D0,                      5.D0,
     4                  7.D0,      9.D0,
     5                       5.D0, 1.D0, 5.D0,
     6                            -1.D0,      5.D0,
     7                                  11.D0,
     8                                        5.D0/

!---------------------------------------------------------------------------
! Allocate storage for the solver ?par and the initial solution vector
!---------------------------------------------------------------------------
      INTEGER length
      PARAMETER (length=128)
      INTEGER ipar(length)
      DOUBLE PRECISION dpar(length),TMP(N,4)
!---------------------------------------------------------------------------
! Some additional variables to use with the RCI (P)CG solver
!---------------------------------------------------------------------------
      CHARACTER MATDES(3)
      DOUBLE PRECISION  solution(N)
      DOUBLE PRECISION expected_sol(N)
      DATA expected_sol/1.D0, 0.D0, 1.D0, 0.D0, 1.D0, 0.D0, 1.D0, 0.D0/
      DOUBLE PRECISION  ONE
      DATA   ONE/1.D0/
      DOUBLE PRECISION DNRM2, Euclidean_norm, temp(N)
      EXTERNAL DNRM2
!---------------------------------------------------------------------------
! Initialize the right hand side through matrix-vector product
!---------------------------------------------------------------------------
       CALL MKL_DCSRSYMV('U', N, A, IA, JA, expected_sol, rhs)
!---------------------------------------------------------------------------
! Initialize the initial guess
!---------------------------------------------------------------------------
       DO I=1, N
         solution(I)=0.D0
       ENDDO
       MATDES(1)='D'
       MATDES(2)='L'
       MATDES(3)='N'
!---------------------------------------------------------------------------
! Initialize the solver
!---------------------------------------------------------------------------
      CALL dcg_init(N, solution,rhs, RCI_request,ipar,dpar,TMP)
      IF (RCI_request .NE. 0 ) GOTO 999
!---------------------------------------------------------------------------
! Set the desired parameters:
! INTEGER parameters:
! set the maximal number of iterations to 100
! LOGICAL parameters:
! run the Preconditioned version of RCI (P)CG with preconditioner C_inverse
! DOUBLE PRECISION parameters
! -
!---------------------------------------------------------------------------
      ipar(5)=100
      ipar(11)=1
!---------------------------------------------------------------------------
! Check the correctness and consistency of the newly set parameters
!---------------------------------------------------------------------------
      CALL dcg_check(N,solution,rhs,RCI_request,ipar,dpar,TMP)
      IF (RCI_request .NE. 0 ) GOTO 999
!---------------------------------------------------------------------------
! Compute the solution by RCI (P)CG solver
! Reverse Communications starts here
!---------------------------------------------------------------------------
1     CALL dcg(N,solution,rhs,RCI_request,ipar,dpar,TMP)
!---------------------------------------------------------------------------
! If RCI_request=0, then the solution was found according to the requested
! stopping tests. In this case, this means that it was found after 100
! iterations.
!---------------------------------------------------------------------------
      IF (RCI_request .EQ. 0) THEN
          GOTO 700
!---------------------------------------------------------------------------
! If RCI_request=1, then compute the vector A*TMP(:,1)
! and put the result in vector TMP(:,2)
!---------------------------------------------------------------------------
      ELSEIF (RCI_request .EQ. 1) THEN
        CALL MKL_DCSRSYMV('U', N, A, IA, JA, TMP, TMP(1,2))
        GOTO 1
!---------------------------------------------------------------------------
! If RCI_request=2, then do the user-defined stopping test: compute the
! Euclidean norm of the actual residual using MKL routines and check if
! it is less than 1.D-8
!---------------------------------------------------------------------------
      ELSEIF (RCI_request .EQ. 2) THEN
        CALL MKL_DCSRSYMV('U', N, A, IA, JA, solution, temp)
        CALL DAXPY(N,-1.D0,rhs,1,temp,1)
        Euclidean_norm = DNRM2(N,temp,1)
        IF (Euclidean_norm .GT. 1.D-8) THEN
!---------------------------------------------------------------------------
! The solution has not been found yet according to the user-defined stopping
! test. Continue RCI (P)CG iterations.
!---------------------------------------------------------------------------
          GOTO 1
        ELSE
!---------------------------------------------------------------------------
! The solution has been found according to the user-defined stopping test
!---------------------------------------------------------------------------
          GOTO 700
        END IF
!---------------------------------------------------------------------------
! If RCI_request=3, then compute apply the preconditioner matrix C_inverse
! on vector TMP(:,3) and put the result in vector TMP(:,4)
!---------------------------------------------------------------------------
      ELSEIF (RCI_request .EQ. 3) THEN
        CALL MKL_DCSRSV('N', N, ONE, MATDES,
     &             A, JA, IA, IA(2), TMP(1,3), TMP(1,4))
        GOTO 1
      ELSE
!---------------------------------------------------------------------------
! If RCI_request=anything else, then dcg subroutine failed
! to compute the solution vector: solution(N)
!---------------------------------------------------------------------------
        GOTO 999
      ENDIF
!---------------------------------------------------------------------------
! Reverse Communication ends here
! Get the current iteration number
!---------------------------------------------------------------------------
700   CALL dcg_get(N,solution,rhs,RCI_request,ipar,dpar,TMP,
     &             itercount)
!---------------------------------------------------------------------------
! Print solution vector: solution(N) and number of iterations: itercount
!---------------------------------------------------------------------------
      WRITE(*, *) ' The system has been solved '
      WRITE(*, *) ' The following solution obtained '
      WRITE(*,800) (solution(i),i =1,N)
      WRITE(*, *) ' expected solution '
      WRITE(*,800)(expected_sol(i),i =1,N)
800   FORMAT(4(F10.3))
      WRITE(*,900)(itercount)
900   FORMAT(' Number of iterations: ',1(I2))
      DO I=1,N
         expected_sol(I)=expected_sol(I)-solution(I)
      ENDDO

      Euclidean_norm=DNRM2(N,expected_sol,1)

!---------------------------------------------------------------------------
! Release internal MKL memory that might be used for computations
! NOTE: It is important to call the routine below to avoid memory leaks
! unless you disable MKL Memory Manager
!---------------------------------------------------------------------------
      CALL MKL_FREEBUFFERS

      IF (itercount.EQ.expected_itercount .AND.
     1                                   Euclidean_norm.LE.1.0D-12) THEN
         WRITE( *,'(A,A)') 'This example has successfully PASSED',
     1 ' through all steps of computation!'
         STOP 0
      ELSE
         WRITE( *,'(A,A,A,I5,A,A,A,E12.5,A)') 'This example may have',
     1 ' FAILED as either the number of iterations differs from the',
     2 ' expected number of iterations ',expected_itercount,' or the',
     3 ' computed solution differs much from the expected solution (',
     4 'Euclidean norm is ',Euclidean_norm,'), or both.'
         STOP 1
      ENDIF
!---------------------------------------------------------------------------
! Release internal MKL memory that might be used for computations
! NOTE: It is important to call the routine below to avoid memory leaks
! unless you disable MKL Memory Manager
!---------------------------------------------------------------------------
999   WRITE( *,'(A,A)') 'This example FAILED as the solver has',
     1 ' returned the ERROR code', RCI_request
      CALL MKL_FREEBUFFERS
      STOP 1

      END

Please find attached a

Juan_Diego_E_ — Fri, 06 Feb 2015 23:54:21 GMT

Please find attached a picture on which you can see what I have added to my project and the error that I have when trying to compile it.

The file name is cg_jacobi

mecej4 — Sat, 07 Feb 2015 03:31:20 GMT

The file name is cg_jacobi_precon.f (if you are using the Windows explorer, you may want to set the folder options to make file suffixes visible). The file suffix is .f, as it should be for a fixed format Fortran source file. There are two errors in the listing that you posted: change the two instances of MKL_FREEBUFFERS to MKL_FREE_BUFFERS.

The screenshots that you attached are not useful because only a part of the source code is shown. At any rate, the source code is provided by Intel in the MKL examples directory or in the examples_core_f.zip file and you can use the provided file without modification.

Good afternoon,

Juan_Diego_E_ — Sat, 07 Feb 2015 19:16:43 GMT

Good afternoon,

Thank you very much mecej4. I will try and change those things you mentioned and check if it works now. In my case I would need to solve a system of equations a lot of times one after the other (it is a transients problem), so everytime I will need to call the function for: initialize, check, and solve? I guess this calls consume a lot of time dont they?

Thank you.

Let us represent the

mecej4 — Sat, 07 Feb 2015 19:59:40 GMT

Let us represent the equations by A.x = b. From what you said, x and b vary with time. But how about A? If A is constant, you only need to factorize it once and use the solution phase of Pardiso from then on, as many time steps as you wish.

If A varies, as well, then the question is whether only certain parts of A change and if you know the exact nature of that change. If the values in A change but the locations of the nonzero elements do not change with time, then all the instances of A(t) are structurally the same, and some economy can be achieved by knowing that.

The more you know the nature of your problem, the better will you be able to choose an appropriate solution strategy.

Good afternoon,

Juan_Diego_E_ — Sat, 07 Feb 2015 20:08:28 GMT

Good afternoon,

Yes, exactly. In my case both A changes (because the main diagonal terms has the time variable there) and B changes as well because this will depend on the previous time-step calculation. The structure of the matrix is exactly the same. So in this case which calls should I do to Pardiso (if I decicde to use this solver and no RCI ISS) at each time step?

I actually just began another forum because I am also trying to figure out which is the best solver for my simulations. Can you please see the forum: PARDISO vs LAPACK Performance ? If you see the code there, you can notice that I am implementing Pardiso with phase 13 (analyze, factorize and solve) at each time step. I dont know if maybe this is the reason why Pardiso is taking almost twice the time vs LAPACK.

Thank you very much mecej4 for your reply. I would appreciate if you can let me know what is the best implementation of Pardiso in this case (please look at the other forum, or if you want I can just post everything here as well)

Regards,

J.D

Pardiso takes a 'phase'

mecej4 — Sat, 07 Feb 2015 20:33:50 GMT

Pardiso takes a 'phase' argument which is a two digit number. The first digit is the starting phase and the second digit is the terminating phase. Thus, 13 is equivalent to 11, 22, 33 in sequence. Of these three, 22 is the most time consuming. For the problem that you described, phase 11 needs to be done only once. The other two phases have to be performed each time step. The savings from doing phase 11 only once may not be large, but definitely worthwhile since the cost of reorganizing your code to do so is very little. In your code in the other thread, you have use 11, 22 and then 13, which amounts to doing 11 and 22 twice each step -- unnecessary and inefficient. If your matrix remains structurally constant, there is no need to deallocate Pardiso memory between time steps.

As to the choice between Lapack and Pardiso, there is no contest. Most Lapack linear solvers are for dense matrices. If your matrix is derived by discretizing a partial differential equation, it will be sparse, possibly banded. If your matrix is banded, Lapack does have some suitable routines, but again Pardiso is probably going to be faster.

There other sparse solvers available (SuperLU, WSMP, UMFPack, etc.) and it is worthwhile for you to learn about their characteristics and ascertain which sparse solver best suits your application.

Good morning mecej4,

Juan_Diego_E_ — Mon, 09 Feb 2015 14:36:49 GMT

Good morning mecej4,

Thank you very much for your observations. I will modify the code accordingly and check if it actually ran faster with Pardiso. If it is not the case, do you think one of the reasons might be that Lapack is actually better for matrices not that big? maybe if the size of my matrices are on the order of 20000 to 45000 I would expect Pardiso to be faster?

Thank you.

Regards,

J.D.

Good afternoon,

Juan_Diego_E_ — Tue, 10 Feb 2015 01:07:01 GMT

Good afternoon,

I am now having trouble to link pardiso in platform 64 bit.

Can anybody please help me? Please find attached a snap shot of how I have everything set up.

Please look at the warning that I get (I dont know why, and how to get rid of it) and the severe error 157.

Thank you very much,

Regards,

J.D.

if (ParTLapF) then
                        allocate(pt(64), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate pt(64), in boundaries1 module ***"
                        do i = 1,64
                           pt( i )%DUMMY =  0 
                        end do
                        maxfct = 1
                        mnum = 1
                        mtype = 2
                        
                        allocate(nzele((size(Tz))+(size(Tth))+(size(Tr))+(Nz*Nth*Nr)), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate nzele((size(Tz))+(size(Tth))+(size(Tr))+(Nz*Nth*Nr)), in boundaries1 module ***"
                        allocate(colnzele((size(Tz))+(size(Tth))+(size(Tr))+(Nz*Nth*Nr)), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate colnzele((size(Tz))+(size(Tth))+(size(Tr))+(Nz*Nth*Nr)), in boundaries1 module ***"
                        allocate(innzele((Nz*Nr*Nth)+1), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate innzele((Nz*Nr*Nth)+1), in boundaries1 module ***"
                        cta = 0
                        ctt = 0
                        !$OMP PARALLEL DO 
                        do i = 1,(Nz*Nr*Nth)
                            do j = 1,(Nz*Nr*Nth)
                                if ((j >= i) .AND. (abs(TMnn(i,j)) > 0.0)) then
                                    ctt = ctt + 1
                                    nzele(ctt) = TMnn(i,j)
                                    colnzele(ctt) = j
                                    if (j == i) then
                                        cta = cta + 1
                                        innzele(cta) = ctt
                                    end if
                                end if
                            end do
                        end do
                        !$OMP END PARALLEL DO 
                        innzele(size(innzele)) = (size(Tz))+(size(Tth))+(size(Tr))+(Nz*Nth*Nr)+1
                                           
                        allocate(perm(Nz*Nr*Nth), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate perm(Nz*Nr*Nth), in boundaries1 module ***"
                        do i = 1,(Nz*Nr*Nth)
                            perm(i) = 1 
                        end do
                        
                        nrhs = 1
                        
                        do i = 1,64
                           iparm(i) = 0
                        end do 

                        iparm(1) = 1 ! no solver default
                        iparm(2) = 2 ! fill-in reordering from METIS
                        iparm(4) = 0 ! no iterative-direct algorithm 0 or 62 (iterative preconditioning)
                        iparm(5) = 0 ! no user fill-in reducing permutation
                        iparm(6) = 1 ! =0 solution on the first n compoments of x
                        iparm(8) = 0 ! numbers of iterative refinement steps
                        iparm(10) = 8 ! perturbe the pivot elements with 1E-13
                        iparm(11) = 0 ! use nonsymmetric permutation and scaling MPS
                        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
                        iparm(14) = 0 ! Output: number of perturbed pivots
                        iparm(24) = 0 ! Parallel factorization control: 1 -> two level factorization algorithm. 0 -> classic algorithm for factorization
                        iparm(25) = 0 ! 0 -> parallel algorithm for the solve step. 1 -> sequential forward and backward solve
                                       
                        msglvl = 0  !not printing statistical information
                        allocate(XXX(Nz*Nr*Nth), STAT = AllocateStatus)
                        if (AllocateStatus /= 0) STOP "*** Not enough memory to allocate XXX(Nz*Nr*Nth), in boundaries1 module ***"
                        
                        E = 0
                        
                    end if
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                    
                            phase = 13 
                            call pardiso_D_64 (pt, maxfct, mnum, mtype, phase, Nz*Nr*Nth, nzele, innzele, colnzele, perm, nrhs, iparm, msglvl, RHSnn, XXX, E)

Hi,

Alexander_K_Intel2 — Tue, 10 Feb 2015 01:15:04 GMT

Hi,

Could you please replace pardiso_D_64 to pardiso and set iparm(27) and msglvl to 1? There could be a problem in matrix structure

Thanks,

Alex

Hi Alex,

Juan_Diego_E_ — Tue, 10 Feb 2015 02:09:08 GMT

Hi Alex,

Thanks for your prompt reply.

This is what I got when I modify what you mentioned. What should I do?

J.D.

Hi Alex,

Juan_Diego_E_ — Tue, 10 Feb 2015 02:31:26 GMT

Hi Alex,

I decided to run it again with pardiso_D_64, but setting the other two parameters you mentioned to 1. And this is the error that I got. What does it mean? how can I fix it?

Thank you very much.

J.D.

Hi,

Alexander_K_Intel2 — Tue, 10 Feb 2015 02:52:12 GMT

Hi,

Which name of include file you use in your example?

About this problem - matrix checker show the issue in your ia array, namely ia(33930) is equal to zero. Pardiso decide that size of your matrix is greater than 33928 and its values is incorrect. What value of n you set in pardiso?

Thanks,

Alex

Hi Alex,

Juan_Diego_E_ — Thu, 12 Feb 2015 00:40:15 GMT

Hi Alex,

Thank you for your reply. I checked that and I corrected the mistake, which actually was because I implemented a OMP PARALLEL DO loop, and it did something I did not want to, so once I performed the a normal do loop, I got the answer.

Thank you and regards,

J.D.

Good afternoon,

Juan_Diego_E_ — Thu, 12 Feb 2015 00:47:18 GMT

Good afternoon,

I would like to ask about the FGMRES with ILU0 Preconditioner solver. I have a code on which I implement LAPACK, PARDISO, and now I would like to check the performance with this iterative solver. LAPACK and PARDISO gave me exactly the same results (as expected) and Pardiso performed 15% to 20% faster than LAPACK for large systems (34000 elements) and around 50-100 time steps. However, when I implemented the FGMRES solver, it actually gave me a completely wrong solution. So I went one step back and first made the example to work, and it did. When I tried to implement the Pardiso example into FGMRES, it actually again gave me a completely different solution. I would like to know, why is this happening? Do I have something wrong? In addition I would like to know what are the best parameters (IPAR and DPAR) I should use for a ill-posed symmetric large (35000 or more) system when using FGMRES?

Please find attached the corresponding documents (Pardiso example and FGMRES example).

Thank you for your reply.

Regards,

J.D.

You attached only the .SLN

mecej4 — Thu, 12 Feb 2015 14:45:31 GMT

You attached only the .SLN files, which are of no use -- we need the source code, along with any data and a list of the compiler options used. You talk of "the Pardiso example", but the MKL examples/solverf/source directory contains six Pardiso examples. Similarly, there are three GMRES related examples. Which of these did you modify?

With all these uncertainties, and without seeing the code, it is not possible to respond.

Good morning mecej4,

Juan_Diego_E_ — Thu, 12 Feb 2015 15:02:51 GMT

Good morning mecej4,

Perfect. Please find attached the source codes for both examples (I am implementing modules, as I was trying to emulate the situation in which I really need to implement this solvers). The first two source codes correspond to PardisoEX (PardisoEX.f90 and input.f90) and the last two source codes correspond to FGMRESex3 (FGMRESex3.f90 and input_0.f90).

I am comparing: pardiso_sym_f90 and dcsrilu0_exampl2

I am using 64 bit platform release mode. Machine: 32 GB RAM, Intel Core i7 3.4 GHz. 4 cores = 8 threads. Visual Studio 2008. Fortran 90. Intel Compiler 11.1-072

On both cases - Properties:

* I included directory in Fortran/general: "C:\Program Files (x86)\Intel\Compiler\11.1\072\mkl\include"

* Fortran/general/Optimization: Maximize speed

* Fortran/Preprocessor/Default Include and Use Path: Source File directory

* Fortran/Preprocessor/OpenMP Conditional Compilation: Yes

* Fortran/Libraries/run time library: Multithread

* Fortran/Libraries/Use Intel Math Kernel Library: Parallel

*Linker/General/Link library dependencies: Yes

* Linker/Manifest File/Generate manifest file: No

In the Pardiso example I also manually included the following:

* In the Header Files Folder: mkl_lapack95.lib

* In the Source Files Folder: lapack.f90 and mkl_pardiso.f90 in addition to the two source files attached to this message.

In the FGMRES example I also manually included the following:

* In the Header Files Folder: mkl_lapack95.lib and mkl_lapack95_lp64.lib

* In the Source Files Folder: lapack.f90, mkl.fi and mkl_rci.fi

Thank you very much for the support.

Regards,

J.D

Here is one possible source

mecej4 — Thu, 12 Feb 2015 19:46:29 GMT

Here is one possible source of error: you apparently are comparing solutions for linear equations with a symmetric indefinite matrix. In the GMRES version of your code, you initialize the (rowidx, col, val) arrays for the upper triangular part of the matrix. The solver routines that you call, however, expect the lower triangular part to be filled in, as well, because these routines (with names containing GE rather than SYM) have to work with unsymmetric matrices too. In effect, then, the Pardiso version gives the solution for the symmetric matrix whereas the GMRES version gives the solution for an upper triangular, unsymmetric matrix. Naturally, the two solutions will not agree.

Good afternoon mecej4,

Juan_Diego_E_ — Thu, 12 Feb 2015 20:18:41 GMT

Good afternoon mecej4,

Thank you very much for your input! I did not think about that. I will modify the inputs accordingly and will come back to you with the results. As you said, I am actually trying to solve a system of equations that comes from a finite difference discretization of a PDE, which results in a symmetric diagonally dominant matrix (with high condition number as the CG iterative solver from MKL did not work giving me the error that the matrix was almost indefinite).

Thanks again.

Regards,

J.D

Good morning mecej4,

Juan_Diego_E_ — Sat, 14 Feb 2015 17:00:04 GMT

Good morning mecej4,

I was busy working these days in other things, and right now I tried what you told me with the Intel examples, and it actually worked! That was my mistake, I was inputting the triangular part of the matrix to the FMGRES solver whereas it should have been the full symmetric matrix.

Thanks a lot for your help!

Regards,

J.D