2Algebra for linear equation system with bordered band matrix.
10# Bordered band matrix.
12# \author Claus Kleinwort, DESY, 2021 (Claus.Kleinwort@desy.de)
74 nSizeBand = nSize - nBorder
86 self.
__mixed = np.zeros((nBorder, nSizeBand))
88 self.
__band = np.zeros((nBand + 1, nSizeBand))
100 for i
in range(len(aIndex)):
101 iIndex = aIndex[i] - 1
102 for j
in range(i + 1):
103 jIndex = aIndex[j] - 1
104 if (iIndex < nBorder):
105 self.
__border[iIndex, jIndex] += aMatrix[i, j]
106 if (iIndex != jIndex):
107 self.
__border[jIndex, iIndex] += aMatrix[i, j]
108 elif (jIndex < nBorder):
109 self.
__mixed[jIndex, iIndex - nBorder] += aMatrix[i, j]
111 nBand = iIndex - jIndex
112 self.
__band[nBand, jIndex - nBorder] += aMatrix[i, j]
128 aMatrix = np.empty((nSize, nSize))
129 for i
in range(nSize):
130 iIndex = aIndex[i] - 1
131 for j
in range(i + 1):
132 jIndex = aIndex[j] - 1
133 iMax = max(iIndex, jIndex)
134 iMin = min(iIndex, jIndex)
136 aMatrix[i, j] = self.
__border[iIndex, jIndex]
137 elif (iMin < nBorder):
138 aMatrix[i, j] = self.
__mixed[iMin, iMax - nBorder]
142 aMatrix[i, j] = self.
__band[iBand, iMin - nBorder]
145 aMatrix[j, i] = aMatrix[i, j]
151 print(
" block part ")
153 for i
in range(nRow):
155 print(
" mixed part ")
156 for i
in range(nRow):
157 print(
" row ", i, self.
__mixed[i])
160 for i
in range(nRow):
161 print(
" diagonal ", i, self.
__band[i])
180 nRow = self.__numBand + 1
182 auxVec = np.copy(self.__band[0]) * 16.0
183 for i
in range(nCol):
184 if ((self.__band[0, i] + auxVec[i]) != self.__band[0, i]):
185 self.__band[0, i] = 1.0 / self.__band[0, i]
187 self.__band[0, i] = 0.0
188 for j
in range(min(nRow, nCol - i) - 1):
189 rxw = self.__band[j + 1, i] * self.__band[0, i]
190 for k
in range(min(nRow, nCol - i) - j - 1):
191 self.__band[k, i + j + 1] -= self.__band[k + j + 1, i] * rxw
192 self.__band[j + 1, i] = rxw
200 nRow = self.__numBand + 1
202 aSolution = np.copy(aRightHandSide)
203 for i
in range(nCol):
204 for j
in range(min(nRow, nCol - i) - 1):
205 aSolution[j + i + 1] -= self.__band[j + 1, i] * aSolution[i]
206 for i
in range(nCol - 1, -1, -1):
207 rxw = self.__band[0, i] * aSolution[i]
208 for j
in range(min(nRow, nCol - i) - 1):
209 rxw -= self.__band[j + 1, i] * aSolution[j + i + 1]
218 nRow = self.__numBand + 1
220 inverseBand = np.zeros((nRow, nCol))
221 for i
in range(nCol - 1, -1, -1):
222 rxw = self.__band[0, i]
223 for j
in range(i, max(0, i - nRow + 1) - 1, -1):
224 for k
in range(j + 1, min(nCol, j + nRow)):
225 rxw -= inverseBand[abs(i - k), min(i, k)] * self.__band[k - j, j]
226 inverseBand[i - j, j] = rxw
238 nBand = self.__numBand
239 aBand = np.empty((nBand + 1, nCol))
240 for i
in range(nCol):
241 for j
in range(max(0, i - nBand), i + 1):
242 aBand[i - j, j] = np.dot(anArray[i], np.dot(aSymArray, anArray[j]))
246 nBorder = self.__numBorder
249 aSolution = np.empty(nBorder + nCol)
252 if ((self.__band[0] <= 0.0).any()):
253 raise ZeroDivisionError(
"Band matrix not positive definite")
256 auxMat = np.empty((nBorder, nCol))
258 for i
in range(nBorder):
262 auxVec = aRightHandSide[:nBorder] - np.dot(auxMat, aRightHandSide[nBorder:])
263 invBorder = np.linalg.inv(self.__border - np.dot(self.__mixed, auxMatT))
264 aSolution[:nBorder] = np.dot(invBorder, auxVec)
266 aSolution[nBorder:] =
solveBand(aRightHandSide[nBorder:]) \
267 -np.dot(auxMatT, aSolution[:nBorder])
269 self.__border = invBorder
270 self.__mixed = np.dot(-invBorder, auxMat)
(Symmetric) Bordered Band Matrix.
__numBand
size of border; int
def addBlockMatrix(self, aIndex, aMatrix)
Add (compressed) block to BBmatrix:
__band
band part C; matrix(float)
__numBorder
size of border; int
__border
border part B; matrix(float)
def getBlockMatrix(self, aIndex)
Retrieve (compressed) block from BBmatrix:
def printMatrix(self)
Print BBmatrix.
__mixed
mixed part M; matrix(float)
def __init__(self, nSize, nBorder=1, nBand=7)
Create new BBmatrix.
__numSize
size of matrix; int
def solveAndInvertBorderedBand(self, aRightHandSide)
Solve linear equation A*x=b system with BBmatrix A, calculate BB part of inverse of A.
__numCol
size of band part of matrix; int
def solveBand(aRightHandSide)
Solve linear equation system for band part.
def bandOfAVAT(anArray, aSymArray)
Calculate band part of A*V*A^T.
def invertBand()
Invert band part.