INCLUDE 'mkl_pardiso.f90'                                                                                                                                                                                      
PROGRAM pardiso_sym_f90                                                                                                                                                                                        
USE mkl_pardiso
 implicit none

 real,dimension(:),allocatable :: sparse_values,b,x
!double precision,dimension(:),allocatable :: sparse_values,b,x
 integer,dimension(:),allocatable :: sparse_columns, sparse_rowIndex, perm

 type(mkl_pardiso_handle) :: pt(64)
 integer :: maxfct,mnum,mtype,phase,n,msglvl,ierr,nrhs
 integer :: iparm(64)

 integer :: i,ij,iteration

 do i=1,64
   pt(i)%DUMMY =  0
 enddo

 iparm = 0
 
 maxfct= 1
  mnum = 1
 mtype = 11
     n = 27
 msglvl= 1
  nrhs = 1
 
 iparm(1) = 1
 iparm(2) = 2
 iparm(3) = 1
 iparm(4) = 0
 iparm(5) = 0
 iparm(6) = 0
 iparm(7) = 0
 iparm(8) = 2
 iparm(9) = 0
 iparm(10) = 13
 iparm(11) = 1
 iparm(12) = 0
 iparm(13) = 1
 iparm(14) = 0
 iparm(15) = 0
 iparm(16) = 0
 iparm(17) = 0
 iparm(18) = -1
 iparm(19) = -1
 iparm(20) = 0
!iparm(27) = 1    ! matrix checker
 iparm(28) = 1    ! single precision

!write(nout,*) "IPARM:", iparm

 allocate(perm(n)); perm=0
 allocate(sparse_rowIndex(n+1)); sparse_rowindex=0
 allocate(b(n),x(n)); b=0.0; x=0.0
 allocate(sparse_values(135),sparse_columns(135)); sparse_values=0.0; sparse_columns=0.0

 b = (/ 2.6884846E+17, 2.0396429E+17,  4.2347395E+15, 1.8505861E+16,&
        1.4292538E+16, 9.6922534E+14,  2.4841678E+13, 3.3282958E+13, 2.5818447E+13,&
        5.3769434E+17, 4.2185506E+17,  7.6176066E+15, 3.7177104E+16, 2.7938754E+16,&
        1.6816193E+15, 9.0916562E+12,  1.3865087E+13, 4.1425329E+13, 2.3148379E+17,&
        1.7876580E+17, 3.9469260E+15,  1.6763004E+16, 1.3027155E+16, 9.2449486E+14,&
        2.2259052E+13, 3.0544388E+13,  2.7990157E+13 /)

 sparse_rowIndex = (/1,  5, 10, 14, 19,&
                    25, 30, 34, 39, 43, 48,&
                    54, 59, 65, 72, 78, 83,&
                    89, 94, 98,103,107,112,&
                   118,123,127,132,136 /)

 sparse_columns = (/      1,  2,  4, 10,  1,&
                      2,  3,  5, 11,  2,  3,&
                      6, 12,  1,  4,  5,  7,&
                     13,  2,  4,  5,  6,  8,&
                     14,  3,  5,  6,  9, 15,&
                      4,  7,  8, 16,  5,  7,&
                      8,  9, 17,  6,  8,  9,&
                     18,  1, 10, 11, 13, 19,&
                      2, 10, 11, 12, 14, 20,&
                      3, 11, 12, 15, 21,  4,&
                     10, 13, 14, 16, 22,  5,&
                     11, 13, 14, 15, 17, 23,&
                      6, 12, 14, 15, 18, 24,&
                      7, 13, 16, 17, 25,  8,&
                     14, 16, 17, 18, 26,  9,&
                     15, 17, 18, 27, 10, 19,&
                     20, 22, 11, 19, 20, 21,&
                     23, 12, 20, 21, 24, 13,&
                     19, 22, 23, 25, 14, 20,&
                     22, 23, 24, 26, 15, 21,&
                     23, 24, 27, 16, 22, 25,&
                     26, 17, 23, 25, 26, 27,&
                     18, 24, 26, 27/)



 do iteration=1,3
   if(iteration == 1) sparse_values = (/ 1.8199702E+17, -1.0985973E+16, -8.1912757E+14, -8.6409330E+15,&
                        -1.1592385E+16,  1.4588944E+17, -1.6598948E+14, -5.6100655E+14, -4.7898050E+15,&
                        -2.0017022E+14,  3.1136505E+15, -1.9205374E+13, -5.9795118E+13, -1.3057714E+15,&
                         1.4310660E+16, -4.9571076E+14, -2.7139611E+11, -3.1083547E+14, -8.6506768E+14,&
                        -5.1366198E+14,  1.1123487E+16, -2.9914416E+13, -3.2975526E+11, -2.3459568E+14,&
                        -2.7025094E+13, -3.1866762E+13,  7.2196575E+14, -1.8579852E+11, -1.1541145E+13,&
                        -3.2787746E+11,  1.8697716E+13, -1.6858925E+12, -3.9060259E+11, -3.8703976E+11,&
                        -1.6760384E+12,  2.5464555E+13, -8.2631020E+11, -7.2688127E+11, -2.3000133E+11,&
                        -5.9658456E+11,  1.8608376E+13, -2.8582235E+11, -8.2937740E+15,  3.5702905E+17,&
                        -2.6536052E+16, -2.0436672E+15, -6.9282499E+15, -4.5801451E+15, -2.7880778E+16,&
                         3.0043138E+17, -3.1712260E+14, -1.2732341E+15, -3.8033412E+15, -5.6776276E+13,&
                        -3.7229968E+14,  5.6106350E+15, -3.3883381E+13, -4.8613879E+13, -3.0381478E+14,&
                        -3.1016079E+15,  2.9225968E+16, -1.0906780E+15, -9.5587705E+10, -2.5677595E+14,&
                        -2.2928706E+14, -1.9219051E+15, -1.1272249E+15,  2.2035849E+16, -5.5743823E+13,&
                        -1.3566280E+11, -1.9474403E+14, -1.1203748E+13, -4.7868887E+13, -5.8888569E+13,&
                         1.2579326E+15, -3.0559918E+11, -9.7358394E+12, -4.0293109E+11, -1.1471707E+11,&
                         8.5987846E+12, -2.0592403E+12, -3.6577552E+11, -7.3020709E+11, -1.5613501E+11,&
                        -2.0062270E+12,  1.3477062E+13, -1.4252345E+12, -6.7511668E+11, -2.7899835E+11,&
                        -3.7727915E+11, -8.9305330E+11,  2.9937778E+13, -2.5662416E+11, -7.3386785E+15,&
                         1.5517277E+17, -9.9426045E+15, -7.5088987E+14, -4.0517559E+15, -1.0494409E+16,&
                         1.2732889E+17, -1.5332513E+14, -5.1321598E+14, -5.1979443E+13, -1.8264756E+14,&
                         2.8940416E+15, -1.8055436E+13, -2.6498432E+14, -1.1933087E+15,  1.2956211E+16,&
                        -4.5460179E+14, -2.1820277E+11, -2.0089857E+14, -7.8911716E+14, -4.7076781E+14,&
                         1.0125458E+16, -2.8321281E+13, -2.6898506E+11, -1.0108888E+13, -2.5336910E+13,&
                        -2.9912136E+13,  6.8683701E+14, -1.9152165E+11, -3.6779567E+11, -2.5767510E+11,&
                         1.6941785E+13, -1.7538799E+12, -6.8959640E+11, -3.1031525E+11, -1.7352727E+12,&
                         2.3639041E+13, -9.2134723E+11, -2.6378419E+11, -2.3387315E+11, -6.5677925E+11,&
                         2.0127326E+13 /)

   if(iteration == 2) sparse_values = (/ 3.0497395E+17, -2.5613430E+16, -2.4096675E+15, -1.0264903E+16,&
                        -1.3849243E+16,  2.4234588E+17, -4.1531039E+14, -1.5931377E+15, -4.3839220E+15,&
                        -1.7329567E+14,  4.8868895E+15, -5.5966930E+13, -4.1442550E+13, -9.7703070E+14,&
                         2.2914033E+16, -1.1035711E+15, -7.7642124E+11, -2.6614256E+14, -6.2921533E+14,&
                        -5.3159361E+14,  1.7885936E+16, -7.4839755E+13, -9.3599898E+11, -1.9779440E+14,&
                        -1.6530695E+13, -2.1230944E+13,  1.1362630E+15, -7.7001110E+11, -7.1963928E+12,&
                        -1.5059164E+11,  2.9406548E+13, -3.5750541E+12, -6.7386809E+11, -1.6881000E+11,&
                        -2.4075690E+12,  4.1018041E+13, -1.9080616E+12, -1.2092768E+12, -1.5144532E+10,&
                        -8.2258048E+11,  2.9104496E+13, -1.8615850E+11, -2.3114554E+16,  6.0231226E+17,&
                        -6.1906456E+16, -5.7562182E+15, -5.5354956E+15, -1.2851870E+16, -3.6267460E+16,&
                         5.0192132E+17, -7.9318154E+14, -3.4733092E+15, -2.1649934E+15, -1.6049888E+14,&
                        -3.8159768E+14,  8.7908322E+15, -9.9741061E+13, -1.7294369E+13, -8.4024653E+14,&
                        -2.5390241E+15,  4.6810372E+16, -2.3679567E+15, -1.7452489E+11, -1.2373749E+14,&
                        -6.2520644E+14, -1.3842096E+15, -1.2149284E+15,  3.5301455E+16, -1.3986758E+14,&
                        -2.6452035E+11, -9.2442158E+13, -3.2230691E+13, -2.7889351E+13, -4.4319491E+13,&
                         1.9744561E+15, -1.2603995E+12, -2.6164171E+12, -1.0294930E+12, -2.3219693E+11,&
                         1.4565046E+13, -4.3841282E+12, -4.7363876E+11, -1.8840568E+12, -2.7574734E+11,&
                        -3.5770500E+12,  2.3590588E+13, -3.1848585E+12, -8.2325615E+11, -7.9882145E+11,&
                        -3.5376316E+10, -1.5946056E+12,  4.6660843E+13, -7.7362913E+10, -1.8507794E+16,&
                         2.5553293E+17, -2.3659199E+16, -2.2484199E+15, -1.0488552E+16, -1.3097261E+16,&
                         2.0915280E+17, -3.8735666E+14, -1.4803319E+15, -1.3491955E+14, -1.7552201E+14,&
                         4.5022545E+15, -5.3581974E+13, -6.6648235E+14, -7.7086315E+14,  2.0515626E+16,&
                        -1.0190081E+15, -5.7569850E+11, -4.9738486E+14, -4.9069783E+14, -4.9760998E+14,&
                         1.6104289E+16, -7.1406188E+13, -7.1203016E+11, -2.6820221E+13, -1.3504847E+13,&
                        -2.0892810E+13,  1.0712847E+15, -7.8354992E+11, -8.3301637E+11, -1.9072367E+11,&
                         2.6595125E+13, -3.7064296E+12, -1.6087163E+12, -2.0509627E+11, -2.5140407E+12,&
                         3.7952227E+13, -2.1463368E+12, -6.9892820E+11, -3.1654500E+10, -9.2923272E+11,&
                         3.1211781E+13/)

   if(iteration == 3) sparse_values = (/ 3.1020985E+17, -2.4146940E+16, -2.1968032E+15, -1.0630049E+16,&
                        -1.4169095E+16,  2.2792357E+17, -3.9137316E+14, -1.4124705E+15, -4.2461123E+15,&
                        -2.5345015E+14,  4.6284779E+15, -5.1097062E+13, -4.0183219E+13, -9.6843386E+14,&
                         2.1218966E+16, -9.7970278E+14, -7.7146096E+11, -2.5265583E+14, -6.0402454E+14,&
                        -5.5577625E+14,  1.5947305E+16, -6.6088637E+13, -8.8817710E+11, -1.8314715E+14,&
                        -1.6257566E+13, -3.4204870E+13,  1.0365083E+15, -7.0787052E+11, -6.9038679E+12,&
                        -1.8506283E+11,  3.0249444E+13, -3.6607025E+12, -8.1247935E+11, -2.1401515E+11,&
                        -2.8427986E+12,  4.0046829E+13, -1.8596013E+12, -1.3092627E+12, -1.8674792E+10,&
                        -9.4394509E+11,  2.7508076E+13, -1.8443366E+11, -2.3771729E+16,  6.3843459E+17,&
                        -6.0378852E+16, -5.3123017E+15, -5.9758409E+15, -1.2098926E+16, -3.7585413E+16,&
                         4.8635701E+17, -7.7988922E+14, -3.1131311E+15, -2.2053970E+15, -1.5142300E+14,&
                        -5.2297777E+14,  8.5245190E+15, -9.2258297E+13, -1.7459978E+13, -7.7167167E+14,&
                        -2.5212207E+15,  4.3910479E+16, -2.1396915E+15, -2.0338260E+11, -1.2375781E+14,&
                        -5.5201258E+14, -1.3444418E+15, -1.2692116E+15,  3.1872486E+16, -1.2595582E+14,&
                        -2.9939794E+11, -9.1343578E+13, -2.9288806E+13, -2.7773192E+13, -6.7652374E+13,&
                         1.8287424E+15, -1.1925013E+12, -2.7521088E+12, -1.2716220E+12, -3.0258489E+11,&
                         1.6516554E+13, -5.0188832E+12, -5.5676056E+11, -1.9836806E+12, -3.5691227E+11,&
                        -4.6337475E+12,  2.4912287E+13, -3.2328951E+12, -9.1093657E+11, -7.5204375E+11,&
                        -3.6856676E+10, -1.8503116E+12,  4.5168795E+13, -7.9658402E+10, -1.9826735E+16,&
                         2.8068455E+17, -2.3839661E+16, -2.1102090E+15, -1.0183662E+16, -1.4106058E+16,&
                         2.0915494E+17, -3.8676744E+14, -1.3491130E+15, -1.3040333E+14, -2.5076093E+14,&
                         4.4525413E+15, -5.0248332E+13, -6.2180423E+14, -7.8546805E+14,  1.9657415E+16,&
                        -9.3759753E+14, -6.0587455E+11, -4.4623777E+14, -4.8732175E+14, -5.3075370E+14,&
                         1.4850881E+16, -6.5598134E+13, -7.1351199E+11, -2.4770874E+13, -1.3571612E+13,&
                        -3.2975480E+13,  1.0101377E+15, -7.5603522E+11, -1.0671994E+12, -2.1115614E+11,&
                         2.8289268E+13, -3.9466851E+12, -1.7914308E+12, -2.3478239E+11, -3.1616399E+12,&
                         3.8514771E+13, -2.1625491E+12, -6.7463289E+11, -3.0868685E+10, -1.1446616E+12,&
                         3.0518366E+13/)

   if(iteration == 1) then

     phase=11

     call pardiso(pt, maxfct, mnum, mtype, phase, n, sparse_values, sparse_rowIndex, &
                  sparse_columns, perm, nrhs, iparm, msglvl, b, x, ierr)

     write(*,*) 'Reordering completed ... '
     if (ierr .ne. 0) then
       write(*,*) 'The following ERROR was detected: ', ierr
       stop 1
     endif
     write(*,*) 'Number of nonzeros in factors = ',iparm(18)
     write(*,*) 'Number of factorization MFLOPS = ',iparm(19)
     phase=23
     call pardiso(pt, maxfct, mnum, mtype, phase, n, sparse_values, sparse_rowIndex, &
                  sparse_columns, perm, nrhs, iparm, msglvl, b, x, ierr)

   else
     phase=23
     iparm(4)=31
     call pardiso(pt, maxfct, mnum, mtype, phase, n, sparse_values, sparse_rowIndex, &
                  sparse_columns, perm, nrhs, iparm, msglvl, b, x, ierr)
   endif

   write(*,'(a)') "Solve completed ..."
   write(*,'(a)') "The solution of the system is: "
   do ij = 1,n
     write(*,'(a,i0,a,f9.6)') " x [",ij,"] = ",x(ij)
   enddo

 enddo

END PROGRAM
