[Pkg-electronics-commits] [gnucap] 13/47: matrix: fbsubt and s

[felix salfelder]

This is an automated email from the git hooks/post-receive script.

felix-guest pushed a commit to branch master
in repository gnucap.

commit 9a1ffcbbf9cb9338a28fa993bcfef73e9cfd86cb
Author: Felix Salfelder <felix at salfelder.org>
Date:   Wed Sep 14 00:59:40 2016 +0100

    matrix: fbsubt and s
---
 include/m_matrix.h | 56 +++++++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 43 insertions(+), 13 deletions(-)

diff --git a/include/m_matrix.h b/include/m_matrix.h
index a18519e..6cb5140 100644
--- a/include/m_matrix.h
+++ b/include/m_matrix.h
@@ -150,6 +150,7 @@ public:
   int		size()const		{return _size;}
   double 	density();
   T 	d(int r, int  )const	{return *(_diaptr[r]);}
+  T const& s(int r, int c)const;
 private:
   T 	u(int r, int c)const	{return _colptr[c][r];}
   T 	l(int r, int c)const	{return _rowptr[r][-c];}
@@ -157,7 +158,6 @@ private:
   T&	u(int r, int c);
   T&	l(int r, int c);
   T&	m(int r, int c);
-  //T&	s(int r, int c);
 public:
   void		load_diagonal_point(int i, T value);
   void		load_point(int i, int j, T value);
@@ -169,6 +169,7 @@ public:
   void		lu_decomp();
   void		fbsub(T* v) const;
   void		fbsub(T* x, const T* b, T* c = NULL) const;
+  void		fbsubt(T* v) const;
 };
 /*--------------------------------------------------------------------------*/
 // private implementations
@@ -445,10 +446,10 @@ T& BSMATRIX<T>::m(int r, int c)
   return (c>=r) ? u(r,c) : l(r,c);
 }
 /*--------------------------------------------------------------------------*/
-/* s: general matrix entry access.
+/* s: general matrix entry access (read-only)
  * It is known that the location is strictly in bounds,
  *   but it is not known whether the location actually exists.
- * If access is attempted to a non-allocated location, 
+ * If access is attempted to a non-allocated location,
  *   it returns a reference to a shared zero variable.
  *   Writing to this zero is not prohibited,
  *   but will corrupt the matrix in a known and testable way.
@@ -457,9 +458,8 @@ T& BSMATRIX<T>::m(int r, int c)
  *   Writing to trash is allowed and encouraged,
  *   but reading it gives a number not useful for anything.
  */
-#if 0
 template <class T>
-T& BSMATRIX<T>::s(int row, int col)
+T const& BSMATRIX<T>::s(int row, int col)const
 {untested();
   assert(_lownode);
   assert(0 <= col);
@@ -469,28 +469,27 @@ T& BSMATRIX<T>::s(int row, int col)
   assert(_zero == 0.);
 
   if (col == row) {untested();
-    return d(row,col);
-  }else if (col > row) {untested();	/* above the diagonal */
+    return d(row, col);
+  }else if (col > row) {untested();    /* above the diagonal */
     if (row == 0) {untested();
       return _trash;
     }else if (row < _lownode[col]) {untested();
       return _zero;
     }else{untested();
-      return u(row,col);
+      return u(row, col);
     }
-  }else{untested();			/* below the diagonal */
+  }else{untested();                    /* below the diagonal */
     assert(col < row);
     if (col == 0) {untested();
       return _trash;
     }else if (col < _lownode[row]) {untested();
       return _zero;
     }else{untested();
-      return l(row,col);
+      return l(row, col);
     }
   }
   unreachable();
 }
-#endif
 /*--------------------------------------------------------------------------*/
 template <class T>
 void BSMATRIX<T>::load_point(int i, int j, T value)
@@ -715,7 +714,7 @@ void BSMATRIX<T>::fbsub(T* x, const T* b, T* c) const
     }
 
     int first_nz = ii;
-    for (   ; ii <= size(); ++ii) {		/* forward substitution */
+    for (   ; ii <= size(); ++ii) {            /* forward substitution */
       int low_node = std::max(_lownode[ii], first_nz);
       c[ii] = b[ii];
       for (int jj = low_node; jj < ii; ++jj) {
@@ -727,7 +726,7 @@ void BSMATRIX<T>::fbsub(T* x, const T* b, T* c) const
 
   notstd::copy_n(c, size()+1, x);
 
-  for (int jj = size(); jj > 1; --jj) {		/* back substitution    */
+  for (int jj = size(); jj > 1; --jj) {         /* back substitution    */
     for (int ii = _lownode[jj]; ii < jj; ++ii) {
       x[ii] -= u(ii,jj) * x[jj];
     }
@@ -737,6 +736,37 @@ void BSMATRIX<T>::fbsub(T* x, const T* b, T* c) const
   // x[0]==0 eliminates a lot of "if" statements
 }
 /*--------------------------------------------------------------------------*/
+/* fbsubt: forward and back substitution with implicitly transposed matrix Ut Lt x = v
+ * v = right side vector, changed in place to solution vector
+ * GS:
+ * this method s used to solve system A_t X = B (_t - transposed)
+ * which corresponds to adjoint system
+ * (LU)_t then transforms to U_t L_t
+ *  added: Gennadiy Serdyuk <gserdyuk at gserdyuk.com>
+ */
+template <class T>
+void BSMATRIX<T>::fbsubt(T* v) const
+{ untested();
+  assert(_lownode);
+  assert(v);
+
+  // forward substitution
+  for (unsigned ii = 1; ii <= size(); ++ii) { untested();
+    for (unsigned jj = _lownode[ii]; jj < ii; ++jj) { untested();
+      v[ii] -= u(jj,ii) * v [jj];
+    }
+  }
+
+  // back substitution
+  for (unsigned jj = size(); jj > 1; --jj) { untested();
+    v[jj] /= d(jj,jj);
+    for (unsigned ii = _lownode[jj]; ii < jj; ++ii) { untested();
+      v[ii] -= l(jj,ii) * v[jj];
+    }
+  }
+  v[1]/=d(1,1);
+}
+/*--------------------------------------------------------------------------*/
 /*--------------------------------------------------------------------------*/
 #endif
 

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/pkg-electronics/gnucap.git