GeneralBrokenLines V03-01-02
gblpy3
gblnum.py
Go to the documentation of this file.
1'''
2Algebra for linear equation system with bordered band matrix.
3
4Created on Jul 27, 2011
5
6@author: kleinwrt
7'''
8
9## \file
10# Bordered band matrix.
11#
12# \author Claus Kleinwort, DESY, 2021 (Claus.Kleinwort@desy.de)
13#
14# \copyright
15# Copyright (c) 2011 - 2023 Deutsches Elektronen-Synchroton,
16# Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
17# This library is free software; you can redistribute it and/or modify
18# it under the terms of the GNU Library General Public License as
19# published by the Free Software Foundation; either version 2 of the
20# License, or (at your option) any later version. \n\n
21# This library is distributed in the hope that it will be useful,
22# but WITHOUT ANY WARRANTY; without even the implied warranty of
23# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24# GNU Library General Public License for more details. \n\n
25# You should have received a copy of the GNU Library General Public
26# License along with this program (see the file COPYING.LIB for more
27# details); if not, write to the Free Software Foundation, Inc.,
28# 675 Mass Ave, Cambridge, MA 02139, USA.
29
30import numpy as np
31
32
33
65class BorderedBandMatrix(object):
66
67
73 def __init__(self, nSize, nBorder=1, nBand=7):
74 nSizeBand = nSize - nBorder
75
76 self.__numSize = nSize
77
78 self.__numBorder = nBorder
79
80 self.__numBand = 0 # actual band width
81
82 self.__numCol = nSizeBand
83
84 self.__border = np.zeros((nBorder, nBorder))
85
86 self.__mixed = np.zeros((nBorder, nSizeBand))
87
88 self.__band = np.zeros((nBand + 1, nSizeBand))
89# print " new BBM ", self.__border__.shape, self.__mixed__.shape, self.__band__.shape
90
91
98 def addBlockMatrix(self, aIndex, aMatrix):
99 nBorder = self.__numBorder
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]
110 else:
111 nBand = iIndex - jIndex
112 self.__band[nBand, jIndex - nBorder] += aMatrix[i, j]
113
114 self.__numBand = max(self.__numBand, nBand)
115 return self.__numBand
116
117
124 def getBlockMatrix(self, aIndex):
125 nBorder = self.__numBorder
126 nBand = self.__numBand
127 nSize = len(aIndex)
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)
135 if (iMax < nBorder):
136 aMatrix[i, j] = self.__border[iIndex, jIndex]
137 elif (iMin < nBorder):
138 aMatrix[i, j] = self.__mixed[iMin, iMax - nBorder]
139 else:
140 iBand = iMax - iMin
141 if iBand <= nBand:
142 aMatrix[i, j] = self.__band[iBand, iMin - nBorder]
143 else:
144 aMatrix[i, j] = 0.
145 aMatrix[j, i] = aMatrix[i, j]
146 return aMatrix
147
148
150 def printMatrix(self):
151 print(" block part ")
152 nRow = self.__numBorder
153 for i in range(nRow):
154 print(" row ", i, self.__border[i])
155 print(" mixed part ")
156 for i in range(nRow):
157 print(" row ", i, self.__mixed[i])
158 nRow = self.__numBand + 1
159 print(" band part ")
160 for i in range(nRow):
161 print(" diagonal ", i, self.__band[i])
162
163
170 def solveAndInvertBorderedBand(self, aRightHandSide):
171
172#============================================================================
173
180 nRow = self.__numBand + 1
181 nCol = self.__numCol
182 auxVec = np.copy(self.__band[0]) * 16.0 # save diagonal elements
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]
186 else:
187 self.__band[0, i] = 0.0 # singular
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
193
194
199 def solveBand(aRightHandSide):
200 nRow = self.__numBand + 1
201 nCol = self.__numCol
202 aSolution = np.copy(aRightHandSide)
203 for i in range(nCol): # forward substitution
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): # backward substitution
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]
210 aSolution[i] = rxw
211 return aSolution
212
213
218 nRow = self.__numBand + 1
219 nCol = self.__numCol
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
227 rxw = 0.0
228 return inverseBand
229
230
236 def bandOfAVAT(anArray, aSymArray):
237 nCol = self.__numCol
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]))
243 return aBand
244
245# setup
246 nBorder = self.__numBorder
247# nRow = self.__numBand +1
248 nCol = self.__numCol
249 aSolution = np.empty(nBorder + nCol)
250# decompose band
252 if ((self.__band[0] <= 0.0).any()):
253 raise ZeroDivisionError("Band matrix not positive definite")
254
255 if (nBorder > 0):
256 auxMat = np.empty((nBorder, nCol))
257# solve for mixed part
258 for i in range(nBorder):
259 auxMat[i] = solveBand(self.__mixed[i])
260 auxMatT = auxMat.T
261# solve for border
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)
265# solve for band part
266 aSolution[nBorder:] = solveBand(aRightHandSide[nBorder:]) \
267 -np.dot(auxMatT, aSolution[:nBorder])
268# parts of inverse
269 self.__border = invBorder
270 self.__mixed = np.dot(-invBorder, auxMat)
271 self.__band = invertBand() + bandOfAVAT(auxMatT, invBorder)
272 else:
273# solve for band part (only)
274 aSolution = solveBand(aRightHandSide)
275 self.__band = invertBand()
276 return aSolution
(Symmetric) Bordered Band Matrix.
Definition: gblnum.py:65
__numBand
size of border; int
Definition: gblnum.py:80
def addBlockMatrix(self, aIndex, aMatrix)
Add (compressed) block to BBmatrix:
Definition: gblnum.py:98
__band
band part C; matrix(float)
Definition: gblnum.py:88
__numBorder
size of border; int
Definition: gblnum.py:78
__border
border part B; matrix(float)
Definition: gblnum.py:84
def getBlockMatrix(self, aIndex)
Retrieve (compressed) block from BBmatrix:
Definition: gblnum.py:124
def printMatrix(self)
Print BBmatrix.
Definition: gblnum.py:150
__mixed
mixed part M; matrix(float)
Definition: gblnum.py:86
def __init__(self, nSize, nBorder=1, nBand=7)
Create new BBmatrix.
Definition: gblnum.py:73
__numSize
size of matrix; int
Definition: gblnum.py:76
def solveAndInvertBorderedBand(self, aRightHandSide)
Solve linear equation A*x=b system with BBmatrix A, calculate BB part of inverse of A.
Definition: gblnum.py:170
__numCol
size of band part of matrix; int
Definition: gblnum.py:82
def solveBand(aRightHandSide)
Solve linear equation system for band part.
Definition: gblnum.py:199
def decomposeBand()
Definition: gblnum.py:179
def bandOfAVAT(anArray, aSymArray)
Calculate band part of A*V*A^T.
Definition: gblnum.py:236
def invertBand()
Invert band part.
Definition: gblnum.py:217