Methods for solving Ax = b

Author: rbotafogo

lib/blis/base/level-1v.rb

@@ -89,6 +89,8 @@ def self.copyv(vx, vy)
   end
 
   #------------------------------------------------------------------------------------
+  # rho := conjx(x)^T * conjy(y)
+  #
   # Performs a dotv operation between two vectors.  No error checking is done. Vectors
   # can be row or column vectors that no error will be issued.
   # @param vecx [MDArray] an MDArray configured as a vector: dimension [x, 1] or [1, x]

lib/linalg/algebra.rb

@@ -57,14 +57,17 @@ def self.forward_solve(mLU, vb)
 
     def self.back_substitution(mLU, vb)
 
-      mUpart = mLU.part_by(:six, row_dir: :rl, column_dir: :bt, filter: 0b001000)
+      mUpart = mLU.part_by(:six, row_dir: :rl, column_dir: :bt, filter: 0b100010)
       vbpart = vb.part_by(:row, column_dir: :bt, filter: 0b11)
 
       loop do
-        beta1, bottom = vbpart.next
-        beta1.pp
-        bottom.pp
-        
+        beta1, b2 = vbpart.next
+        v11, u12t = mUpart.next
+        #  b1     1    beta1    -1      u12t         b2
+        # rho := beta * rho + alpha * conjx(x)^T * conjy(y)
+        # dotxv(alpha, vx, vy, beta, rho)
+        Blis.dotxv(-1, u12t, b2, 1, beta1)
+        beta1[0, 0] /= v11[0, 0]
       end
 
     end

lib/linalg/decomposition.rb

@@ -31,15 +31,16 @@ class LinAlg
     # @param a [MDArray]
     #--------------------------------------------------------------------------------------
 
-    def self.lu(a)
+    def self.lu(mA)
 
-      apart = a.part_by(:six, row_dir: :lr, column_dir: :tb, filter: 0b101011)
+      apart = mA.part_by(:six, row_dir: :lr, column_dir: :tb, filter: 0b101011)
 
       loop do
         diag, l21, a12, a22 = apart.next
         Blis.scalv(1.0/diag[0, 0], l21)
         Blis.ger(-1, a22, l21, Blis.trnsp(a12))
       end
+      mA
 
     end
 

test/test_decomposition.rb

@@ -63,9 +63,15 @@ class MDArrayTest < Test::Unit::TestCase
     #--------------------------------------------------------------------------------------
 
     should "LU decompose a matrix" do
-=begin
+
+      lu_a = MDArray.double([4, 4],
+                            [2, 0, 1, 2,
+                             -1, -1, 2, 1,
+                             2, 1, 1, -1,
+                             -2, -1, 1, -2])
+      
       MDArray::LinAlg.lu(@a)
-      @a.pp
+      assert_equal(true, lu_a.identical(@a))
 
       MDArray::LinAlg.lu(@b)
       @b.pp
@@ -84,7 +90,6 @@ class MDArrayTest < Test::Unit::TestCase
       b = MDArray.double([3, 1], [2, 10, -2])
       MDArray::LinAlg.forward_solve(@d, b)
       b.pp
-=end
 
       MDArray::LinAlg.lu(@c)
       @c.pp

test/test_solvers.rb

@@ -0,0 +1,97 @@
+# -*- coding: utf-8 -*-
+
+##########################################################################################
+# Copyright © 2017 Rodrigo Botafogo. All Rights Reserved. Permission to use, copy, modify, 
+# and distribute this software and its documentation, without fee and without a signed 
+# licensing agreement, is hereby granted, provided that the above copyright notice, this 
+# paragraph and the following two paragraphs appear in all copies, modifications, and 
+# distributions.
+#
+# IN NO EVENT SHALL RODRIGO BOTAFOGO BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, 
+# INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF 
+# THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF RODRIGO BOTAFOGO HAS BEEN ADVISED OF THE 
+# POSSIBILITY OF SUCH DAMAGE.
+#
+# RODRIGO BOTAFOGO SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
+# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE 
+# SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS PROVIDED "AS IS". 
+# RODRIGO BOTAFOGO HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, 
+# OR MODIFICATIONS.
+##########################################################################################
+
+require "test/unit"
+require 'shoulda'
+
+require '../config' if @platform == nil
+require 'mdarray-blis'
+
+
+class MDArrayTest < Test::Unit::TestCase
+
+  context "LinAlg" do
+
+    #--------------------------------------------------------------------------------------
+    #
+    #--------------------------------------------------------------------------------------
+
+    setup do
+      
+      @a = MDArray.double([4, 4],
+                          [2, 0, 1, 2,
+                           -2, -1, 1, -1,
+                           4, -1, 5, 4,
+                           -4, 1, -3, -8])
+
+      @b = MDArray.double([4, 1], [2, 10, -2, 5])
+
+      
+      @c = MDArray.double([4, 4],
+                          [2, 0, 1, 2,
+                           -2, -1, 1, -1,
+                           4, -1, 5, 4,
+                           -4, 1, -3, -8])
+
+      @b_orig = MDArray.double([4, 1], [2, 10, -2, 5])
+      
+    end
+
+    #--------------------------------------------------------------------------------------
+    #
+    #--------------------------------------------------------------------------------------
+
+    should "forward solve an equation " do
+
+      # perform LU factorization of a
+      MDArray::LinAlg.lu(@a)
+
+      a_b = MDArray.double([4, 1], [2, 12, -18, 39])
+      back_sub = MDArray.double([4, 1], [39.25, -106.5, -37.5, -19.5])
+
+      MDArray::LinAlg.forward_solve(@a, @b)
+      assert_equal(true, a_b.identical(@b))
+
+      MDArray::LinAlg.back_substitution(@a, @b)
+      assert_equal(true, back_sub.identical(@b))
+
+      vy = MDArray.double([4, 1])
+      # def self.gemv_dot(alpha, mA, vx, beta, vy)
+      Blis.gemv_dot(1, @c, @b, 1, vy)
+      assert_equal(true, @b_orig.identical(vy))
+
+    end
+    
+    #--------------------------------------------------------------------------------------
+    #
+    #--------------------------------------------------------------------------------------
+
+    should "solve an equation " do
+
+      MDArray::LinAlg.back_substitution(
+        @a, MDArray::LinAlg.forward_solve(
+          MDArray::LinAlg.lu(@a), @b))
+
+    end
+    
+  end
+
+end