Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/sparse/linalg/_isolve/utils.py: 14%

1__docformat__ = "restructuredtext en" 

2 

3__all__ = [] 

4 

5 

6from numpy import asanyarray, asarray, array, zeros 

7 

8from scipy.sparse.linalg._interface import aslinearoperator, LinearOperator, \ 

9 IdentityOperator 

10 

11_coerce_rules = {('f','f'):'f', ('f','d'):'d', ('f','F'):'F', 

12 ('f','D'):'D', ('d','f'):'d', ('d','d'):'d', 

13 ('d','F'):'D', ('d','D'):'D', ('F','f'):'F', 

14 ('F','d'):'D', ('F','F'):'F', ('F','D'):'D', 

15 ('D','f'):'D', ('D','d'):'D', ('D','F'):'D', 

16 ('D','D'):'D'} 

17 

18 

19def coerce(x,y): 

20 if x not in 'fdFD': 

21 x = 'd' 

22 if y not in 'fdFD': 

23 y = 'd' 

24 return _coerce_rules[x,y] 

25 

26 

27def id(x): 

28 return x 

29 

30 

31def make_system(A, M, x0, b): 

32 """Make a linear system Ax=b 

33 

34 Parameters 

35 ---------- 

36 A : LinearOperator 

37 sparse or dense matrix (or any valid input to aslinearoperator) 

38 M : {LinearOperator, Nones} 

39 preconditioner 

40 sparse or dense matrix (or any valid input to aslinearoperator) 

41 x0 : {array_like, str, None} 

42 initial guess to iterative method. 

43 ``x0 = 'Mb'`` means using the nonzero initial guess ``M @ b``. 

44 Default is `None`, which means using the zero initial guess. 

45 b : array_like 

46 right hand side 

47 

48 Returns 

49 ------- 

50 (A, M, x, b, postprocess) 

51 A : LinearOperator 

52 matrix of the linear system 

53 M : LinearOperator 

54 preconditioner 

55 x : rank 1 ndarray 

56 initial guess 

57 b : rank 1 ndarray 

58 right hand side 

59 postprocess : function 

60 converts the solution vector to the appropriate 

61 type and dimensions (e.g. (N,1) matrix) 

62 

63 """ 

64 A_ = A 

65 A = aslinearoperator(A) 

66 

67 if A.shape[0] != A.shape[1]: 

68 raise ValueError(f'expected square matrix, but got shape={(A.shape,)}') 

69 

70 N = A.shape[0] 

71 

72 b = asanyarray(b) 

73 

74 if not (b.shape == (N,1) or b.shape == (N,)): 

75 raise ValueError(f'shapes of A {A.shape} and b {b.shape} are ' 

76 'incompatible') 

77 

78 if b.dtype.char not in 'fdFD': 

79 b = b.astype('d') # upcast non-FP types to double 

80 

81 def postprocess(x): 

82 return x 

83 

84 if hasattr(A,'dtype'): 

85 xtype = A.dtype.char 

86 else: 

87 xtype = A.matvec(b).dtype.char 

88 xtype = coerce(xtype, b.dtype.char) 

89 

90 b = asarray(b,dtype=xtype) # make b the same type as x 

91 b = b.ravel() 

92 

93 # process preconditioner 

94 if M is None: 

95 if hasattr(A_,'psolve'): 

96 psolve = A_.psolve 

97 else: 

98 psolve = id 

99 if hasattr(A_,'rpsolve'): 

100 rpsolve = A_.rpsolve 

101 else: 

102 rpsolve = id 

103 if psolve is id and rpsolve is id: 

104 M = IdentityOperator(shape=A.shape, dtype=A.dtype) 

105 else: 

106 M = LinearOperator(A.shape, matvec=psolve, rmatvec=rpsolve, 

107 dtype=A.dtype) 

108 else: 

109 M = aslinearoperator(M) 

110 if A.shape != M.shape: 

111 raise ValueError('matrix and preconditioner have different shapes') 

112 

113 # set initial guess 

114 if x0 is None: 

115 x = zeros(N, dtype=xtype) 

116 elif isinstance(x0, str): 

117 if x0 == 'Mb': # use nonzero initial guess ``M @ b`` 

118 bCopy = b.copy() 

119 x = M.matvec(bCopy) 

120 else: 

121 x = array(x0, dtype=xtype) 

122 if not (x.shape == (N, 1) or x.shape == (N,)): 

123 raise ValueError(f'shapes of A {A.shape} and ' 

124 f'x0 {x.shape} are incompatible') 

125 x = x.ravel() 

126 

127 return A, M, x, b, postprocess