Fast solve of linear problems

General
1

Dear GetFEM users,

I would like to fast up as much as possible the execution of a trivial heat conduction problem where the external heat flux is the only changing parameter with time.
In principle, the tangent matrix remains the same whereas the RHS changes at each iteration.

If I add only linear terms they will be assembled only once, which is already a considerable time saving.
On the other hand, I assume that the linsolve procedure is repeated at each iteration, meaning by this the inversion of the tangent matrix, am I right?
For example in case of linsolve_lu, the factorization is repeated each time even thou the tangent matrix is the same?
If so, is there a way (in python) to access the result of the matrix inversion so as to be able to perform just a simple “inverse_matrix*RHS” multiplication at each iteration?

Thank you.

2

Dear Lorenzo,

Your request makes very much sense. The solution depends very much on which solver you want to use. But if we e.g. talk about mumps, currently if you use

dx = gf.linsolve_mumps(md.tangent_matrix(), md.rhs())

many times for different rhs, the linsolve_mumps command will every time do the same LU factorization before solving. It will also use the unsymmetric mumps version which is slower than the symmetric version.

We should redesign this interface to allow something like

fact = gf.factorize_mumps(md.tangent_matrix(), "sym")
dx = gf.linsolve_mumps(fact, md.rhs())
.... # call gf.linsolve_mumps as many times as necessary
delete(fact) # close connection to MUMPS

the difficulty is that we need to create a new python object for “fact” that will store all information related to an active connection to MUMPS.

The factorize_mumps method should look like this in C++

    id.job = JOB_INIT;
    id.par = 1;
    id.sym = sym ? 2 : 0;
    id.comm_fortran = USE_COMM_WORLD;
    mumps_interf<T>::mumps_c(id);

    id.n = int(gmm::mat_nrows(A));
    id.nz = int(AA.irn.size());
    id.irn = &(AA.irn[0]);
    id.jcn = &(AA.jcn[0]);
    id.a = (MUMPS_T*)(&(AA.a[0]));

    id.ICNTL(1) = -1; // output stream for error messages
    id.ICNTL(2) = -1; // output stream for other messages
    id.ICNTL(3) = -1; // output stream for global information
    id.ICNTL(4) = 0;  // verbosity level
//    id.ICNTL(14) += 80; // small boost to the workspace size
    id.ICNTL(31) = 1;   // only factorization, no solution to follow
    id.ICNTL(33) = 1;   // request determinant calculation

    id.job = 4; // analysis (job=1) + factorization (job=2)
    mumps_interf<T>::mumps_c(id);
    mumps_error_check(id);

Then linsolve mumps should do:

    id.rhs = (MUMPS_T*)(&(rhs[0]));
    id.job = 3;
    mumps_interf<T>::mumps_c(id);
    bool ok = mumps_error_check(id);

And when the object is deleted:

    id.job = JOB_END;
    mumps_interf<T>::mumps_c(id);
3

Dear Kostas,
thank you for your very comprehensive reply.

Is this interface redesign something that is already in the implementation plans by the GetFEM development team? If it’s not the case, I could try by myself but my experience with C++ is poor for now and it will not be able to succeed in in the short term.

Thank you.

4

it has been at least in my plans for quite some time. There are 3 python libraries that implement this kind of interface

We should replicate the same logic in the getfem scripting interface. The implementation could be done in the file interface/src/gf_linsolve.cc

But it requires to create a new “mumps context” object, initialized with something like:

  auto mumpsctx = std::make_shared<getfem::some_new_name_for_mumps_context_cpp_object>
                  (arg1, arg2, ...);
  id_type id = store_mumps_ctx_object(shp);
  workspace().set_dependence(id, (const void *)(...));
  out.pop().from_object_id(id, MUMPS_CTX_CLASS_ID);
}

The method store_mumps_ctx_object need to be implemented in interface/src/getfemint.h

and the MUMPS_CTX_CLASS_ID identifier needs to be added to the getfemint_class_id enumeration in interface/src/getfemint.h.

It is not crazy difficult but it requires to understand the logic of the scripting interface in GetFEM, which is not python specific.

5

Thank you Kostas.

I think that for my short term testing purposes I could try out one of the 3 python libraries and input the relevant tangent matrix and RHS vector from GetFEM, just to check the performance improvement for my problem.
I’ve started giving a look to the GetFEM C++ interface to external solvers but I need more time and less work commitments to be able to sort out everything and contribute effectively.

Thank you for your advices.

6

Sounds like a good idea. It should not be difficult to convert the tangent matrix to any format that 3rd party libraries may require. Fore example in the following snippet that I use for eigenmode analysis, getfem.Spmat matrices are converted to scipy.sparse.csc_matrix:

import scipy.sparse.linalg
#import scipy.io

...
II = {}
...
II["u"] = range(*md.interval_of_variable("u").cumsum())
...

    K = gf.Spmat("copy", md.tangent_matrix(), II["u"], II["u"])
    print("K_uu matrix determinant %ge%i" % K.determinant()[0::2])
    K = scipy.sparse.csc_matrix((K.csc_val(),*(K.csc_ind()[::-1])))
    # scipy.io.savemat(f"K_{i}.mat", {'K':K}) # for exporting the matrix to matlab
    # K = scipy.io.loadmat(f"K_{i}.mat")["K"] # for importing the matrix from matlab
    M = gf.asm_generic(mim, 2, "rho(x)*Test2_u.Test_u", -1, md, "select_output", "u", "u")
    M = scipy.sparse.csc_matrix((M.csc_val(),*(M.csc_ind()[::-1])))
    vals,vecs = scipy.sparse.linalg.eigs(K,6,M,sigma=0.)
    vals = vals.real
    vecs = vecs.real.T
    print("EIGENVALUES:", " ".join(map(str,vals)))

    ...
    output = ()
    ...
    for ii,vec in enumerate(vecs):
      output += (mfu, vec, f"eigenmode_{ii}")
    mfout.export_to_vtu("filname.vtu", *output)