topic Alex: in Intel® oneAPI Math Kernel Library

PARDISO and FEAST

JohnNichols — Sun, 06 Mar 2016 05:52:42 GMT

The old 73 Harrison Structural Analysis program written in 66 Fortran has now had the solver taken out and PARDISO put in its place. It was an interesting task as I tried to ensure it was all in MODULES.

MODULE StructureType
    USE SOLVER
    INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
    REAL (KIND=dp) :: g = 9.806, pi = 3.14159265D0

So the interesting question is the use of the standard DP constant, if I use it in each module, they conflict when I call the solver module from the structure type module - it is also useful having the g and pi defined. I can fix it by making DP private in SOLVER, but that means if I use solver elsewhere alone, I have to unprivate DP to declare the main program variables - How am I supposed to do it for real?

      MODULE Solver
      INTEGER, PRIVATE, PARAMETER :: dp = selected_real_kind(15, 307)

      contains

      subroutine linsolve(n, error, nnz, a, ja, ia, b, x)
      implicit none
      include 'mkl_pardiso.fi'

It is fascinating climbing into someone's old code, Harrison created a GK matrix to hold the data for the inversion routine, quite clever, but not easy to unwind when you want to use PARDISO - took a while to sort out his names for the various arrays -- the book gave some clues, but it took a while. Slowly cleaning up the code -- was it worth it -- I think so as he had excellent notes -- and a good sample problem, although he made some mistakes with his units conversions but it means I have a nice front end for PARDISO -- which is somewhat faster than my normal structures package and a simple text file is nice for small structures models.

The FEAST is a real treat, it works now from within the program so I get the solution and the modes, should not be hard to do time stepping, but that can wait. The program takes less than 0.5 seconds to solve a simple problem - my issue is all the writing out from the different codes, Pardios is only 0.015 seconds of time.

STEVE: Can I suggest that if INTEL put PARDISO and FEAST together you would have an excellent routine - not that it is not that hard, but the different methods for storing matrices is a bit of a waste and calling them one after the other is tedious.

do 567 j=1,mn_6
            do 566 L4 = 1,mn_12
                Results(J) = Results(J) + SATT(I,J,L4)*disp(L4)
566 END do
567 end do

The S matrix (location in XYZ) and the A transpose matrix are set up in SAT to do the final beam resultant forces and moments it is quickest to take the SAT results adn store each one in an 3 D array where I is each beam. I is an outside loop, is there a way to simplify this set of loops or am I stuck with them.

One way is to male a module

andrew_4619 — Sun, 06 Mar 2016 09:15:45 GMT

One way is to male a module that just has the kinds and constants in it and then USE that in the other modules.

John, feel free to discuss

Steven_L_Intel1 — Sun, 06 Mar 2016 18:09:59 GMT

John, feel free to discuss your suggestion in the MKL forum.

Quote:andrew_4619 wrote:

JohnNichols — Mon, 07 Mar 2016 02:08:51 GMT

andrew_4619 wrote:

One way is to male a module that just has the kinds and constants in it and then USE that in the other modules.

I was thinking that is the best way - but it means cascaded modules -- life is not perfect - except maybe if you like Lisp

Hi All,

Alexander_K_Intel2 — Mon, 07 Mar 2016 16:18:14 GMT

Hi All,

I am both in Pardiso and Feast but looks like missed something in thread :)

What do you mean under this:

Can I suggest that if INTEL put PARDISO and FEAST together you would have an excellent routine - not that it is not that hard, but the different methods for storing matrices is a bit of a waste and calling them one after the other is tedious.

Thanks,

Alex

Alex:

JohnNichols — Wed, 09 Mar 2016 22:03:40 GMT

Alex:

I have a very old 1973 structural analysis package that I found in Harrison's book.

The package creates a structures stiffness matrix, quite simple for linear structures. Program is 900 lines long.

The solver from 73 is essentially Gaussian Elimination and it is ok. We normally use STRAND7 to solve our structural problems, although there are many such programs, all big kludgy and bloated. The Italian guy was using STRAND to look at the beam of interest to us at the moment.

The first thing we wanted was the modes, so I added FEAST to the end of 73 structures program, and it worked well, although it has taken a while to work out if I set L to the Number of Degrees of Freedom of the matrix I get the exact solution for the max number of modes. You have L as a parameter, which needs to be changed if you want to pass in the right number.

I also wanted to check the results from the '73 program so I pulled out the '73 solver and threw in PARDISO, it was actually quite easy to create the stiffness matrix in the right form for PARDISO. Now I probably have the world's fastest solver tied to the world's fastest Eigen routine and in less than a second I can do what takes a lot longer time with STRAND 7 and I can modify the printout as required -- DRAIN3 on steroids.

The next problem to solve is the mass distribution - not to hard.

If I can run the stiffness array A(N,N) into Pardiso and then into the eigensolver and along the way pass in a B(N) matrix and get the solution X(n) directly for the results for the inversion, and send it onto the EIGEN routine without having to fix A in another form -- bloody winner.

Simple really -- about 400 lines of FORTRAN code without checks.

Just a thought

John