Coverage for kwant/linalg/_mumps.pyx : 74%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Copyright 2011-2013 Kwant authors.
2#
3# This file is part of Kwant. It is subject to the license terms in the file
4# LICENSE.rst found in the top-level directory of this distribution and at
5# https://kwant-project.org/license. A list of Kwant authors can be found in
6# the file AUTHORS.rst at the top-level directory of this distribution and at
7# https://kwant-project.org/authors.
9cimport numpy as np
10import numpy as np
11from . cimport cmumps
12from . import cmumps
13from .fortran_helpers import assert_fortran_matvec, assert_fortran_mat
15int_dtype = cmumps.int_dtype
17# Proxy classes for Python access to the control and info parameters of MUMPS
19cdef class mumps_int_array:
20 cdef cmumps.MUMPS_INT *array
22 def __init__(self):
23 self.array = NULL
25 def __getitem__(self, key):
26 return self.array[key - 1]
28 def __setitem__(self, key, value):
29 self.array[key - 1] = value
32# workaround for the fact that cython cannot pass pointers to an __init__()
33cdef make_mumps_int_array(cmumps.MUMPS_INT *array):
34 wrapper = mumps_int_array()
35 wrapper.array = array
36 return wrapper
39cdef class smumps_real_array:
40 cdef cmumps.SMUMPS_REAL *array
42 def __init__(self):
43 self.array = NULL
45 def __getitem__(self, key):
46 return self.array[key - 1]
48 def __setitem__(self, key, value):
49 self.array[key - 1] = value
52cdef make_smumps_real_array(cmumps.SMUMPS_REAL *array):
53 wrapper = smumps_real_array()
54 wrapper.array = array
55 return wrapper
58cdef class dmumps_real_array:
59 cdef cmumps.DMUMPS_REAL *array
61 def __init__(self):
62 self.array = NULL
64 def __getitem__(self, key):
65 return self.array[key - 1]
67 def __setitem__(self, key, value):
68 self.array[key - 1] = value
71cdef make_dmumps_real_array(cmumps.DMUMPS_REAL *array):
72 wrapper = dmumps_real_array()
73 wrapper.array = array
74 return wrapper
77cdef class cmumps_real_array:
78 cdef cmumps.CMUMPS_REAL *array
80 def __init__(self):
81 self.array = NULL
83 def __getitem__(self, key):
84 return self.array[key - 1]
86 def __setitem__(self, key, value):
87 self.array[key - 1] = value
90cdef make_cmumps_real_array(cmumps.CMUMPS_REAL *array):
91 wrapper = cmumps_real_array()
92 wrapper.array = array
93 return wrapper
96cdef class zmumps_real_array:
97 cdef cmumps.ZMUMPS_REAL *array
99 def __init__(self):
100 self.array = NULL
102 def __getitem__(self, key):
103 return self.array[key - 1]
105 def __setitem__(self, key, value):
106 self.array[key - 1] = value
109cdef make_zmumps_real_array(cmumps.ZMUMPS_REAL *array):
110 wrapper = zmumps_real_array()
111 wrapper.array = array
112 return wrapper
114#############################################################
116cdef class zmumps:
117 cdef cmumps.ZMUMPS_STRUC_C params
119 cdef public mumps_int_array icntl
120 cdef public zmumps_real_array cntl
121 cdef public mumps_int_array info
122 cdef public mumps_int_array infog
123 cdef public zmumps_real_array rinfo
124 cdef public zmumps_real_array rinfog
126 def __init__(self, verbose=False, sym=0):
127 self.params.job = -1
128 self.params.sym = sym
129 self.params.par = 1
130 self.params.comm_fortran = -987654
132 cmumps.zmumps_c(&self.params)
134 self.icntl = make_mumps_int_array(self.params.icntl)
135 self.cntl = make_zmumps_real_array(self.params.cntl)
136 self.info = make_mumps_int_array(self.params.info)
137 self.infog = make_mumps_int_array(self.params.infog)
138 self.rinfo = make_zmumps_real_array(self.params.rinfo)
139 self.rinfog = make_zmumps_real_array(self.params.rinfog)
141 # no diagnostic output (MUMPS is very verbose normally)
142 if not verbose:
143 self.icntl[1] = 0
144 self.icntl[3] = 0
146 def __dealloc__(self):
147 self.params.job = -2
148 cmumps.zmumps_c(&self.params)
150 def call(self):
151 cmumps.zmumps_c(&self.params)
153 def _set_job(self, value):
154 self.params.job = value
156 def _get_job(self):
157 return self.params.job
159 job = property(_get_job, _set_job)
161 @property
162 def sym(self):
163 return self.params.sym
165 def set_assembled_matrix(self,
166 cmumps.MUMPS_INT N,
167 np.ndarray[cmumps.MUMPS_INT, ndim=1] i,
168 np.ndarray[cmumps.MUMPS_INT, ndim=1] j,
169 np.ndarray[np.complex128_t, ndim=1] a):
170 self.params.n = N
171 self.params.nz = a.shape[0]
172 self.params.irn = <cmumps.MUMPS_INT *>i.data
173 self.params.jcn = <cmumps.MUMPS_INT *>j.data
174 self.params.a = <cmumps.ZMUMPS_COMPLEX *>a.data
176 def set_dense_rhs(self, np.ndarray rhs):
178 assert_fortran_matvec(rhs)
179 if rhs.dtype != np.complex128:
180 raise ValueError("numpy array must be of dtype complex128!")
182 if rhs.ndim == 1:
183 self.params.nrhs = 1
184 else:
185 self.params.nrhs = rhs.shape[1]
186 self.params.lrhs = rhs.shape[0]
187 self.params.rhs = <cmumps.ZMUMPS_COMPLEX *>rhs.data
189 def set_sparse_rhs(self,
190 np.ndarray[cmumps.MUMPS_INT, ndim=1] col_ptr,
191 np.ndarray[cmumps.MUMPS_INT, ndim=1] row_ind,
192 np.ndarray[np.complex128_t, ndim=1] data):
194 if row_ind.shape[0] != data.shape[0]:
195 raise ValueError("Number of entries in row index and value "
196 "array differ!")
198 self.params.nz_rhs = data.shape[0]
199 self.params.nrhs = col_ptr.shape[0] - 1
200 self.params.rhs_sparse = <cmumps.ZMUMPS_COMPLEX *>data.data
201 self.params.irhs_sparse = <cmumps.MUMPS_INT *>row_ind.data
202 self.params.irhs_ptr = <cmumps.MUMPS_INT *>col_ptr.data
204 def set_schur(self,
205 np.ndarray[np.complex128_t, ndim=2, mode='c'] schur,
206 np.ndarray[cmumps.MUMPS_INT, ndim=1] schur_vars):
208 if schur.shape[0] != schur.shape[1]:
209 raise ValueError("Schur matrix must be squared!")
210 if schur.shape[0] != schur_vars.shape[0]:
211 raise ValueError("Number of Schur variables must agree "
212 "with Schur complement size!")
214 self.params.size_schur = schur.shape[0]
215 self.params.schur = <cmumps.ZMUMPS_COMPLEX *>schur.data
216 self.params.listvar_schur = <cmumps.MUMPS_INT *>schur_vars.data