mumpsincludes.h File Reference

#include "../../shared/Numerics/types.h"

Go to the source code of this file.

Functions
void	SeqDenseMumpsSolve (IssmDouble *uf, int uf_M, int uf_n, IssmDouble *Kff, int Kff_M, int Kff_N, int Kff_m, IssmDouble *pf, int pf_M, int pf_m, Parameters *parameters)

Function Documentation

SeqDenseMumpsSolve()

void SeqDenseMumpsSolve	(	IssmDouble *	uf,
		int	uf_M,
		int	uf_n,
		IssmDouble *	Kff,
		int	Kff_M,
		int	Kff_N,
		int	Kff_m,
		IssmDouble *	pf,
		int	pf_M,
		int	pf_m,
		Parameters *	parameters
	)

Definition at line 296 of file MumpsSolve.cpp.

                                                                                                                                                                      {/*{{{*/
 
   /*Variables*/
   int        *irn_loc = NULL;
   int        *jcn_loc = NULL;
   IssmDouble *a_loc   = NULL;
   IssmDouble *rhs     = NULL;
 
   /*First, some checks*/
   if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
   if (Kff_M!=Kff_m)_error_("stiffness matrix Kff is not serial");
   if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf");
   if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf");
 
   /*Initialize matrix */
   /*figure out number of non-zero entries: */
   int nnz = 0;
   for(int i=0;i<Kff_M;i++){
      for(int j=0;j<Kff_N;j++){
         if(Kff[i*Kff_N+j]!=0)nnz++;
      }
   }
 
   /*Allocate: */
   if(nnz){
      irn_loc = xNew<int>(nnz);
      jcn_loc = xNew<int>(nnz);
      a_loc   = xNew<IssmDouble>(nnz);
   }
 
   /*Populate the triplets: */
   int count=0;
   for(int i=0;i<Kff_M;i++){
      for(int j=0;j<Kff_N;j++){
         if(Kff[i*Kff_N+j]!=0){
            irn_loc[count] = i+1; //fortran indexing
            jcn_loc[count] = j+1; //fortran indexing
            a_loc[count]   = Kff[i*Kff_N+j];
            count++;
         }
      }
   }
 
   /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
// AD performance is sensitive to calls to ensurecontiguous.
// Providing "t" will cause ensurecontiguous to be called.
#ifdef _HAVE_AD_
   rhs=xNew<IssmDouble>(pf_M,"t");
#else
   rhs=xNew<IssmDouble>(pf_M);
#endif
   xMemCpy<IssmDouble>(rhs,pf,Kff_M);
 
   MumpsSolve(Kff_M,nnz,nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
 
   /*Now scatter from cpu 0 to all other cpus*/
   xMemCpy<IssmDouble>(uf,rhs,Kff_M);
 
   /*Cleanup*/
   xDelete<int>(irn_loc);
   xDelete<int>(jcn_loc);
   xDelete<IssmDouble>(a_loc);
   xDelete<IssmDouble>(rhs);
}/*}}}*/