More refined control of factorize/solve.

- permit out-of-core storage of factors
- option to compute determinant during factorization
- option to perform steps of iterative refinement after solve
- option to solve with the transposed.

src/jmumps.c

@@ -127,7 +127,7 @@ void mumps_factorize(void *jmumps, int n, int nz,
 
 // Solve
 // Solve a linear system using the factorization computed previously.
-void mumps_solve(void* jmumps, int nrhs, double* rhs) {
+void mumps_solve(void* jmumps, int nrhs, double* rhs, int* transposed) {
 
   DMUMPS_STRUC_C* pmumps = (DMUMPS_STRUC_C*)jmumps;
   int taskid, ierr;
@@ -156,6 +156,8 @@ void mumps_solve(void* jmumps, int nrhs, double* rhs) {
     pmumps->nrhs = nrhs;
     pmumps->lrhs = pmumps->n;
     pmumps->rhs  = rhs;  // Will be overwritten with the solution.
+    if (transposed != 0)
+      pmumps->icntl[8] = 0;
   }
 
   pmumps->job  = JOB_SOLVE;

src/jmumps.h

@@ -5,7 +5,7 @@ int   mumps_initialize_mpi(void);
 int   mumps_finalize_mpi(void);
 void* mumps_initialize(int, int*, double*);
 void  mumps_factorize(void*, int, int, double*, int*, int*);
-void  mumps_solve(void*, int, double*);
+void  mumps_solve(void*, int, double*, int*);
 void  mumps_get_info(void*, int*, double*);
 void  mumps_finalize(void*);
 void  mumps_alloc(DMUMPS_STRUC_C**);

src/mumps.jl

@@ -47,6 +47,7 @@ type Mumps
   infog   :: Array{Int32,1}
   rinfog  :: Array{Float64,1}
   nnz     :: Int               # Number of nonzeros in factors.
+  det     :: Int
   err     :: Int
 
   # Constructor.
@@ -72,13 +73,35 @@ type Mumps
     infog = zeros(Int32, 40);
     rinfog = zeros(Float64, 20);
 
-    self = new(id, int32(sym), icntl, cntl, 0, Float64, infog, rinfog, 0, 0);
+    self = new(id, int32(sym), icntl, cntl, 0, Float64, infog, rinfog, 0, 0, 0);
     finalizer(self, finalize);  # Destructor.
     return self;
   end
 end
 
 
+# Convenience constructor.
+# Note that every argument is a keyword.
+function Mumps(;sym :: Int=0,
+               det :: Bool=false,       # Compute determinant.
+               verbose :: Bool=false,   # Output intermediate info.
+               ooc :: Bool=false,       # Store factors out of core.
+               itref :: Int=0,          # Max steps of iterative refinement.
+               cntl :: Array{Float64,1}=default_cntl
+               )
+
+  icntl = default_icntl[:];
+  icntl[33] = det ? 1 : 0;
+  if !verbose
+    icntl[1:4] = 0;
+  end
+  icntl[22] = ooc ? 1 : 0;
+  icntl[10] = itref;
+
+  return Mumps(sym, icntl=icntl, cntl=cntl);
+end
+
+
 function finalize(mumps :: Mumps)
   # Terminate a Mumps instance.
   @mumps_call(:mumps_finalize, Void, (Ptr{Void},), mumps.__id);
@@ -117,6 +140,9 @@ function factorize(mumps :: Mumps, A :: SparseMatrixCSC)
   mumps.n = n;
   mumps.valtype = valtype;
   mumps.nnz = mumps.infog[29];
+  if mumps.icntl[33] == 1
+    mumps.det = mumps.rinfog[12] * 2^(mumps.infog[34]);
+  end
   mumps.err = mumps.infog[1];
   return;
 end
@@ -125,7 +151,7 @@ end
 factorize(mumps :: Mumps, A :: Array{Float64}) = factorize(mumps, sparse(A));
 
 
-function solve(mumps :: Mumps, rhs :: Array{Float64})
+function solve(mumps :: Mumps, rhs :: Array{Float64}; transposed :: Bool=false)
   n = size(rhs, 1);
   if n != mumps.n
     error("rhs has incompatible dimension");
@@ -137,8 +163,8 @@ function solve(mumps :: Mumps, rhs :: Array{Float64})
   nrhs = size(rhs, 2);
   x = rhs[:];
   @mumps_call(:mumps_solve, Void,
-              (Ptr{Void}, Int32, Ptr{Float64}),
-              mumps.__id,  nrhs,            x);
+              (Ptr{Void}, Int32, Ptr{Float64},              Int32),
+              mumps.__id,  nrhs,            x, transposed ? 1 : 0);
 
   @mumps_call(:mumps_get_info, Void,
               (Ptr{Void}, Ptr{Int32}, Ptr{Float64}),