Coverage for /usr/share/miniconda3/envs/dolfin/lib/python3.8/site-packages/block/algebraic/trilinos/Amesos.py: 0%

1from __future__ import division 

2 

3from builtins import str 

4from block.block_base import block_base 

5 

6serial_solvers = ['Klu', 'Lapack', 'Umfpack', 'Taucs', 'Superlu'] 

7parallel_solvers = ['Superludist', 'Mumps', 'Dscpack', 'Pardiso', 'Paraklete'] 

8 

9class AmesosSolver(block_base): 

10 __doc__ = \ 

11 """Trilinos interface to direct solvers. The available solvers depend on 

12 your Trilinos installation, but the default (Klu) is usually available. 

13 

14 Serial solvers: %s 

15 

16 Parallel solvers: %s 

17 

18 Call AmesosSolver.query() to list currently available solvers. 

19 """%(serial_solvers, parallel_solvers) 

20 

21 default_solver = 'Klu' 

22 

23 def __init__(self, A, solver=default_solver, parameters={}): 

24 from PyTrilinos import Epetra, Amesos 

25 from dolfin import info 

26 from time import time 

27 self.A = A # Keep reference 

28 T = time() 

29 self.problem = Epetra.LinearProblem() 

30 self.problem.SetOperator(A.down_cast().mat()) 

31 self.solver = Amesos.Factory().Create(solver, self.problem) 

32 if (self.solver == None): 

33 raise RuntimeError('Failed to create Amesos solver: '+solver) 

34 

35 # This prevents a use-after-free crash for the MPI communicator. It 

36 # enforces that the solver is destroyed before A. 

37 self.solver.A = A 

38 

39 self.solver.SetParameters(parameters) 

40 if self.solver is None: 

41 raise RuntimeError("Unknown solver '%s'"%solver) 

42 err = self.solver.SymbolicFactorization() 

43 if err != 0: 

44 raise RuntimeError("Amesos " + solver + " symbolic factorization failed, err=%d"%err) 

45 err = self.solver.NumericFactorization() 

46 if err != 0: 

47 raise RuntimeError("Amesos " + solver + " numeric factorization failed, err=%d"%err) 

48 info('constructed direct solver (using %s) in %.2f s'%(solver,time()-T)) 

49 

50 @staticmethod 

51 def query(which=None): 

52 """Return list of available solver backends, or True/False if given a solver name""" 

53 from PyTrilinos import Amesos 

54 factory = Amesos.Factory() 

55 if which is None: 

56 avail = [] 

57 for s in serial_solvers + parallel_solvers: 

58 if factory.Query(s): 

59 avail.append(s) 

60 return avail 

61 return factory.Query(which) 

62 

63 def matvec(self, b): 

64 from dolfin import GenericVector 

65 if not isinstance(b, GenericVector): 

66 return NotImplemented() 

67 if self.A.size(0) != len(b): 

68 raise RuntimeError( 

69 'incompatible dimensions for Amesos matvec, %d != %d'%(len(self.b),len(b))) 

70 x = self.A.create_vec(dim=1) 

71 

72 self.problem.SetLHS(x.down_cast().vec()) 

73 self.problem.SetRHS(b.down_cast().vec()) 

74 err = self.solver.Solve() 

75 self.problem.SetLHS(None) 

76 self.problem.SetRHS(None) 

77 

78 if err != 0: 

79 raise RuntimeError("Amesos solve failed, err=%d"%err) 

80 return x 

81 

82 def __str__(self): 

83 return '<%s solver for %s>'%(self.__class__.__name__, str(self.A)) 

84 

85class MumpsSolver(AmesosSolver): 

86 """Specialized Amesos-Mumps solver. The matrix_type argument may be 'SPD', 

87 'symmetric' or 'general'. 

88 

89 NOTE 1: This functionality (matrix_type) is documented, but disabled in 

90 trilinos (see Amesos_Control.h) 

91 

92 NOTE 2: The key that is used in the actual trilinos code is "MatrixProperty" 

93 but the documentation says "MatrixType". 

94 """ 

95 def __init__(self, A, matrix_type='general', parameters={}): 

96 new_parameters = {'MatrixType': matrix_type, 'MatrixProperty': matrix_type} 

97 new_parameters.update(parameters) 

98 super(MumpsSolver, self).__init__(A, solver='Mumps', parameters=new_parameters) 

99