Easy matrix

[Luxan]

Eventually, having this code as a library and *.bi file might be convenient.

[Luxan]

Moving along, three files to illustrate use as a library.
We might revisit the code from other contributors, for inclusion in this library.
There's a chronicle of how we constructed these routines that is always available for
scrutiny.

---------------------------((( 1 )))-----------------------------------

Code: Select all

' easy_matx.bas

' On my computer
' cd Downloads/fb/matrices_x/dodicat/merge_1/lib
'

' ' compile with: fbc -lib easy_matx.bas

'   Matrix math
'
'  1. dim as Matrix a
'  2. redim a.m( x1 to x2 , y1 to y2 )
'
'  3. dim as Matrix z(z1 to z2)  ?
'  4. redim z.m(w1 to w2)  ?
'
'     Matrix(x , y) ~= redim (0 to x-1, 0 to y-1) ?


' These declarations from :
'  http://www.rosettacode.org/wiki/Matrix_multiplication#FreeBASIC
' That content is available under
'  Content is available under GNU Free Documentation License 1.2 unless otherwise noted.
'

' Preliminary results checked using Maxima CAS
'
' also at, especially for vector by matrix calculations.
' https://keisan.casio.com/exec/system/15052033860538
'
' https://matrix.reshish.com/multCalculation.php
'
' https://elsenaju.eu/Calculator/matrix-vector-product.htm
'
'
type Matrix
    dim as double m( any , any )
    declare constructor ( )
    declare constructor ( byval x as uinteger , byval y as uinteger )
end type
 
constructor Matrix ( )
end constructor
 
constructor Matrix ( byval x as uinteger , byval y as uinteger )
    redim this.m( x - 1 , y - 1 )
end constructor


type Vector
dim as double v( any )
    declare constructor ( )
    declare constructor ( byval x as uinteger )
end type

constructor Vector ( )
end constructor

constructor Vector ( byval x as uinteger )
    redim this.v( x - 1 )
end constructor
'
' ====================================================================
'
Public operator * ( byref a as Matrix , byref b as Matrix ) as Matrix
'  Only applicable to square matrices or arrays.
' This routine from :
'  http://www.rosettacode.org/wiki/Matrix_multiplication#FreeBASIC
'
    dim as Matrix ret
    dim as uinteger i, j, k
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                for k = 0 to ubound( b.m , 1 )
                    ret.m( i , j ) += a.m( i , k ) * b.m( k , j )
                next k
            next j
        next i
    end if
    return ret
end operator
'
'
Public operator * ( byref a as Matrix , byref b as Vector ) as Vector
'
'   Vector multiplied by a matrix, maybe correct .
'
'
   dim as Vector ret
   dim as uinteger i, j, k
   if ubound( a.m , 2 ) = ubound( b.v , 1 )  then
        redim ret.v( ubound( a.m , 1 ) )
        for i = 0 to ubound( a.m , 1 )
            
                for k = 0 to ubound( b.v , 1 )
                    ret.v( i ) += a.m( i , k ) * b.v( k  )
                next k
            
        next i
    end if   
   return ret
end operator
'
Public operator - ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                ret.m( i , j ) = a.m( i , j ) - b.m( i , j )
            next j
        next i
    end if
    return ret
end operator
 '
Public operator + ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
               ret.m( i , j ) = a.m( i , j ) + b.m( i , j )
            next j
        next i
    end if
    return ret
end operator
'
Public operator / ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    dim as double x, y, z
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                x = a.m( i , j )
                y = b.m( i , j )
                y = sgn(x)*10^32 ' assume y = 0
             if y <> 0 then z = x/y
                ret.m( i , j ) = z
            next j
        next i
    end if
    return ret
end operator
'
' ......................................................................
'
Public sub prt_z(z() as Matrix, idx as Integer)
'
'     Matrix of Matrix, index and print .
'
dim as integer i,j
'
'print
for i=0 to ubound(z(idx).m,1)
  for j=0 to ubound(z(idx).m,2)  
    print z(idx).m(i,j);" ";
  next j
  print
next i
print
'  
end sub
'
' ......................................................................
'
Public sub prt_m(z as Matrix)
'
'  Print elements from rectangular matrix .
'
dim as integer i,j
'

print
print ubound(z.m,1);" ";ubound(z.m,2)
print
for i=0 to ubound(z.m,1)
  for j=0 to ubound(z.m,2)  
    print z.m(i,j);" ";
  next j
  print
next i
print
'  
end sub
'
' ......................................................................
'
Public sub prt_v(z as Vector)
'
'  Print elements from vector .
'
dim as integer i
'
'print
for i=0 to ubound(z.v,1)
    
    print z.v(i);" ";
  
  
next i
print
print
'  
end sub
'
' ......................................................................
'
Public sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
'     Multiply vector by matrix .
'
dim as integer i,j,ub1, ub2
dim as single x,y,s
    
    ub1=ubound(M1.m,1)
    ub2=ubound(M1.m,2)
'    
    redim V3.v(0 to ub1)
'
        for j=0 to ub1
             s=0
         for i=0 to ub2
             x=V1.v(i)
             y=M1.m(j,i)
             s=s+x*y
         next i
          V3.v(j)=s
      next j
 '
end sub

---------------------------((( 2 )))-----------------------------------

Code: Select all

' easy_matx.bi


#inclib "easy_matx"

'  For some reason, may need to duplicate types, but not constructors,
' in this *.bi

type Matrix
    dim as double m( any , any )
    declare constructor ( )
    declare constructor ( byval x as uinteger , byval y as uinteger )
end type
'
type Vector
dim as double v( any )
    declare constructor ( )
    declare constructor ( byval x as uinteger )
end type

'
'
' ---------------------------- operators -------------------------------
'
declare operator * ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator * ( byref a as Matrix , byref b as Vector ) as Vector

declare operator - ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator + ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator / ( byref a as Matrix , byref b as Matrix ) as Matrix

'   The rest from luxan, sciwiseg@gmail.com
'
'  Probably need to make this content
'  available under GNU Free Documentation License 1.2 
'
'
' ............................. procedures .............................
'
declare sub prt_z(z() as Matrix, idx as Integer)
declare sub prt_m(z as Matrix)
declare sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)

declare sub prt_v(z as Vector)
'

---------------------------((( 3 )))-----------------------------------

Code: Select all

' matx_test.bas


'' compile with: fbc matx_test.bas
#include once "easy_matx.bi"

'
' ----------------------------------------------------------------------
'
'    Matrix of matrices
'
dim as Matrix z(1 to 6)
'
'some garbage matrices for demonstration
dim as Matrix a = Matrix(4 , 2)
a.m(0 , 0) = 1 : a.m(0 , 1) = 0
a.m(1 , 0) = 0 : a.m(1 , 1) = 1
a.m(2 , 0) = 2 : a.m(2 , 1) = 3
a.m(3 , 0) = 0.75 : a.m(3 , 1) = -0.5
'
dim as Matrix M1 = Matrix( 3 , 3 )
M1.m(0 , 0) = 1 : M1.m(0 , 1) = 2 : M1.m(0 , 2) = 3 
M1.m(1 , 0) = 4 : M1.m(1 , 1) = 5 : M1.m(1 , 2) = 6
M1.m(2 , 0) = 7 : M1.m(2 , 1) = 8 : M1.m(2 , 2) = 9
'M1.m(3 , 0) = 3.75 : M1.m(3 , 1) = -2.5 : M1.m(3 , 2) = 1.5
'
dim as Matrix M2 = Matrix( 3 , 3 )
M2.m(0 , 0) = 2 : M2.m(0 , 1) = 4 : M2.m(0 , 2) = 6 
M2.m(1 , 0) = 8 : M2.m(1 , 1) = 10 : M2.m(1 , 2) = 12
M2.m(2 , 0) = 3 : M2.m(2 , 1) = 7 : M2.m(2 , 2) = 9
'M2.m(3 , 0) = 5.0 : M2.m(3 , 1) = -0.5 : M2.m(3 , 2) = 2.5
'
'   dim multiple matrices, or use array of matrices.
'
'

z(1)=a
z(2)=M1
z(3)=M2
'z(4)=M2
'
'print
'print " <<<< M1 "; ubound(M1.m,1);"  ";ubound(M1.m,2)
'print " >>>> z "; ubound(z(1).m,1);"  ";ubound(z(1).m,2)

print
print("  M1  ")
prt_m(M1)
print("  M2  ")
prt_m(M2)

'  
'prt_z(z() , 2)
'prt_z(z() , 3)

dim as Matrix eg = M1 * M2
print " M1.M2 "
prt_m(eg)

z(4) = z(2) * z(3)
print "  z(2).z(3)  " 
prt_z(z() , 4)

eg = M1 + M2
print " M1 + M2 "
prt_m(eg)

eg = M1 - M2

print " M1 - M2 "
prt_m(eg)
print

redim  M1.m( 3 , 2 )
M1.m(0 , 0) = 1 : M1.m(0 , 1) = 2 : M1.m(0 , 2) = 3 
M1.m(1 , 0) = 4 : M1.m(1 , 1) = 5 : M1.m(1 , 2) = 6
M1.m(2 , 0) = 7 : M1.m(2 , 1) = 8 : M1.m(2 , 2) = 9
M1.m(3 , 0) = 7 : M1.m(3 , 1) = 1 : M1.m(3 , 2) = 3
z(5) = M1

print("  M1  ")
prt_m(M1)
'
dim as Vector V1 = Vector(3)
V1.v(0)=1:V1.v(1)=2:V1.v(2)=3
print(" V1 ")
prt_v(V1)
dim as Vector V3 = Vector(3)
Vect_x_Matrix(M1, V1 , V3 )
print(" M1.V1 ")
prt_v(V3)

dim as Vector V4 = M1 * V1
prt_v(V4)
'
'
print " z(2) "
prt_z(z() , 2)
print " z(5) "
prt_z(z() , 5)

end
'
' ====================================================================
'

[fxm]

Hi fxm

As requested, here's there updated code, utilizing your suggestions.
Incorporated into the rest of the code like this, is it alright with you.

OK now. The program no longer aborts when executing 'Vect_x_Matrix()' if compiling with the '-exx' option.

[dodicat]

IMHO
matrix rows 1 to something
matrix columns 1 to something

Fb easily allows this one based approach, element(1,1) is logically the first element of a matrix.
(0,0) is a computer thing because many languages only have zero based arrays, and because of this many computer buffs are hard wired to zero based arrays.
But my schooling took place a long time before home computers, like many others in this forum, so I am not shackled in any way.
But it is a personal choice of course.
Matrix division is normally Matrix multiplied by inverse of another matrix.
Dividing each element of one by each element of another (rdivide or maybe ./) is rarely used, so IMHO the / operator is wasted on this type of division.
The .bi file is a good idea, but only after you have completed your project.

[fxm]

Code: Select all

'  For some reason, may need to duplicate types, but not constructors,
' in this *.bi

The body of [Type ... End Type] is the place of declaration for member data and member procedures (including constructors, destructor, member operators,...), just like the 'Declare ...' lines of non-member procedures. Compiler needs these at both library compile time and user program compile time, hence duplication.
On the other hand, the body of definition of all procedures (including member procedures, constructors,destructor, member operators, ...) are to be compiled only once, therefore not duplicated.

All these declarations and only these can be grouped together in a single file, in order to be included without duplication error when compiling the library, then when compiling the user programs.

[dodicat]

I have added an init function to my methods.

Code: Select all

''#cmdline "-exx" ''optional
#define roundresult(x,N) rtrim(rtrim(left(str((x)+(.5*sgn((x)))/(10^(N))),instr(str((x)+(.5*sgn((x)))/(10^(N))),".")+(N)),"0"),".")
#define Cround(x,N) iif(roundresult(x,N)="","0",roundresult(x,N))
const decplaces=6
type matrix
    as double element(any,any)
    declare constructor
    Declare Constructor(As Long,As Long)
    declare sub init cdecl(as long,as long,...)
    Declare Operator Cast() As String
end type

constructor matrix
end constructor

Constructor matrix (r As Long,c As Long)
Redim element(1 To r,1 To c)
End Constructor

Operator matrix.cast() As String 'for printing
Dim As String ans,comma
#define gap(x) string(2*decplaces-len(x)," ")
For a As Integer=1 To Ubound(this.element,1)
      For b As Integer=1 To Ubound(this.element,2)
           if b=1 then ans+="["
            If b=Ubound(this.element,2) Then comma="" Else comma=","
            Var q=Str(Cround(element(a,b),decplaces))
            ans=ans+q+comma'+gap(q)
            if b=Ubound(this.element,2) then ans+="]"
      Next b
      ans=ans+Chr(10)
Next a
Operator= ans
End Operator

sub matrix.init cdecl(r as long,c as long,...)
    redim this.element(1 to r,1 to c)
    Dim args As Cva_List 
    Cva_Start( args,c )
    for row as long=1 to r
        for col as long=1 to c
        this.element(row,col)= (Cva_Arg(args, Double))
        next
    next
     Cva_End( args ) 
    end sub


dim as matrix m
    m.init(4,3,  1.0,2.0,3.0,_     
                 4.0,5.0,6.0,_ 
                -9.0,-5.2,9.0,_
                 0.0,.4,-7.5)

print m

sleep 

 

Rules
parameter 1 rows
parameter 2 columns

All the rest are strict doubles.
You can enter the doubles in row and column style, or just one after the other
Also in my other methods (previous post) in the inverse
Print "Can't do inverse"
return unit(1)
end if
Not unit(0) which crashes on -exx.

[fxm]

For the 'operator * ( byref a as Matrix , byref b as Matrix ) as Matrix', I do not understand the double condition:
if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
For me, only the first condition is necessary:
if ubound( a.m , 2 ) = ubound( b.m , 1 ) then

With this change, we can now consider a vector as a matrix with only one column ('Dim As Matrix Vect = Matrix(3, 1)' for example), and the above operator still works (not the case without this change).

[Luxan]

Hi dodicat

All these routines might be written in an alternate version
that utilizes the matrix indexing that starts at 1.
However there might be other FreeBasic matrix libraries that use
the 0 indexing method; to be determined at some date.
I doubt you were born before computers existed, or matrix
arithmetic on computers; Fortran then IMSL had these ages ago.

Ypur init function, interesting to see this applied;
is this compatible with the other routines.
What is the command line compile option -exx do.

Hi fxm

Your clarification, about *.bi files, types and constructors is
useful, destructors maybe declared and used in a similar way ?
For me, a constructor is reminiscent of the C command malloc .


Your suggestion, about removal of the extraneous if..then
statement has been noted and will be implemented. The routine
containing these was derived from the matrix multiplication routine,
from Rosetta Code; and was coded in a hurry.

I'm mirroring these routines, as much as possible, through
Maxima CAS; there a vector is generated as a one column matrix
and addressed like this : v[i,1] or v[1,i], depending upon
whether v[] has been transposed.

[fxm]

- A constructor does not allocate memory for the instance of the Type (this was done before by 'Dim' for example), but only "constructs" the object (initialization/construction if necessary of data members, then execution of the user code in the body of the constructor).
- Symmetrically, the destructor does not deallocate memory for the instance of the Type (this will be done at the end), but only "destroys" the object (execution of the user code in the body of the destructor, then destruction of data members if necessary).

[dodicat]

I am going to scrap constructors altogether.
I'll keep the structure (type).
Since OOP is not needed for matrix arithmetic, no sense in having too much.
My test for writing a library(dll) -- run it in C and Pascal.
My init is made into a function.
And operator cast() goes because it returns a string.
-exx is used for error handling.
it will show out of bound arrays, but for actual crashes, if it catches them it will print the error message to the console.
But a crash usually means the console vanishes, so you need an ide/editor which shows the compiler output on a console which stays open.
You can also use onerror goto to catch -exx outputs.
For fun:

Code: Select all

Type matrix 
      Dim As Double element(Any,Any)
End Type

function decplaces() byref as long export
    static as long d=5
    return d
    end function

sub show(M as matrix) export
 #define roundresult(x,N) rtrim(rtrim(left(str((x)+(.5*sgn((x)))/(10^(N))),instr(str((x)+(.5*sgn((x)))/(10^(N))),".")+(N)),"0"),".")
#define Cround(x,N) iif(roundresult(x,N)="","0",roundresult(x,N))   
Dim As String ans,comma
For a As Integer=1 To Ubound(M.element,1)
      For b As Integer=1 To Ubound(M.element,2)
           if b=1 then print"[";
            If b=Ubound(M.element,2) Then comma="" Else comma=","
            Var q=Str(Cround(M.element(a,b),decplaces))
            print q+comma;
            if b=Ubound(M.element,2) then print"]"
      Next b
Next a
End sub


    function initmatrix cdecl(r as long,c as long,...) as matrix export
        dim as matrix m
    redim m.element(1 to r,1 to c)
    Dim args As Cva_List 
    Cva_Start( args,c  )
    for row as long=1 to r
        for col as long=1 to c
        m.element(row,col)= (Cva_Arg(args, Double))
        next
    next
     Cva_End( args )
     return m
    end function
    
    
    Operator *(m1 As matrix,m2 As matrix) As matrix export
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element,1) Then
      Print "Can't do multiply"
      Exit Operator
End If
Dim As matrix ans
Redim ans.element(rows,columns)
Dim rxc As Double
For r As Integer=1 To rows
      For c As Integer=1 To columns
            rxc=0
            For k As Integer = 1 To Ubound(m1.element,2)
                  rxc=rxc+m1.element(r,k)*m2.element(k,c)
            Next k
            ans.element(r,c)=rxc
      Next c
Next r
Operator= ans
End Operator


 dim as matrix A=initmatrix(4,4,  -8.0,9.0,8.5,7.6,_
                                   7.0,2.0,0.0,8.2,_
                                   6.6,-7.3,-2.0,7.6,_
                                   9.0,0.0,5.0,1.0)
         show(A)
         
         print
         
         show((initmatrix(4,4,  -8.0,9.0,8.5,7.6,_     
                                 7.0,2.0,0.0,8.2,_
                                 6.6,-7.3,-2.0,7.6,_
                                9.0,0.0,5.0,1.0)) *  ( initmatrix(4,4,  1/-8.0,1/9.0,8.5,3/7.6,_
                                                                        7.0,2.0,0.0,8.2,_
                                                                        6.6,1/-7.3,-2.0,7.6,_
                                                                        9.0,0.0,1/5.0,1.0)) _
                            )
                            
                            
        decplaces()=12
        
        print
        
       show((initmatrix(4,4,  -8.0,9.0,8.5,7.6,_     
                                      7.0,2.0,0.0,8.2,_
                                      6.6,-7.3,-2.0,7.6,_
                                      9.0,0.0,5.0,1.0)) *  ( initmatrix(4,4,  1/-8.0,1/9.0,8.5,3/7.6,_
                                                                                  7.0,2.0,0.0,8.2,_
                                                                                  6.6,1/-7.3,-2.0,7.6,_
                                                                                  9.0,0.0,1/5.0,1.0)) _
                            ) 
                            
    sleep
    
    
                                      
                                                        
                                                                               
                                                                                       
 

[fxm]

@dodicat,

In your Init procedures (with variadic arguments), what is the use of the third argument ('sz' or 'size') given that the array is sized by the first 2 arguments ('r' and 'c') and is fully initialized regardless of the variadic arguments provided.
You can initialize the variadic list (by 'Cva_Start') on the second argument ('c') and remove 'sz' or 'size'.

[dodicat]

Yes, that works OK.
I thought that
Cva_Start( args,something )
The something had to be the actual size as an original parameter, I did experiment a bit with it, but I must have missed that.
Thanks, I have altered it.
I made a dll, here is the C code using it.

Code: Select all

#include<stdio.h>
#include<stdlib.h>
#define multiply _ZmlR6MATRIXS0_ 
struct matrix {
 double element[100][100];  
};

  struct matrix  INITMATRIX(long r,long c,...); 
  void     SHOW(struct matrix m); 
  struct matrix multiply(struct matrix m1,struct matrix m2);
  long *DECPLACES() ;
 
 
 int main()
 {
 struct	matrix test,b,c;
 test=INITMATRIX(2,2,   1.0,.6,
                          7.7,9.3);
 SHOW(test);
 puts(" ");
 puts(" ");
 
 
            b=   INITMATRIX(4,4,  -8.0,9.0,8.5,7.6,     
                                      7.0,2.0,0.0,8.2,
                                      6.6,-7.3,-2.0,7.6,
                                      9.0,0.0,5.0,1.0);
                                      
             c= INITMATRIX(4,4,  1/-8.0,1/9.0,8.5,3/7.6,
                                    7.0,2.0,0.0,8.2,
                                    6.6,1/-7.3,-2.0,7.6,
                                    9.0,0.0,1/5.0,1.0);                          
                                      
 SHOW(b);
  puts(" ");
 puts(" ");
 SHOW(c);
  puts(" ");
 puts(" ");
 *DECPLACES()=8;
 SHOW(multiply(b,c));
 puts(" ");
 puts(" ");
 *DECPLACES()=12;
  SHOW(multiply(b,c));
system("pause");
 
 } 

results:

Code: Select all

[1,0.6]
[7.7,9.3]


[-8,9,8.5,7.6]
[7,2,0,8.2]
[6.6,-7.3,-2,7.6]
[9,0,5,1]


[-0.125,0.11111,8.5,0.39474]
[7,2,0,8.2]
[6.6,-0.13699,-2,7.6]
[9,0,0.2,1]


[188.5,15.94672755,-83.48,142.84210526]
[86.925,4.77777778,61.14,27.36315789]
[3.275,-13.59269406,61.62,-64.85473684]
[40.875,0.31506849,66.7,42.55263158]


[188.5,15.946727549467,-83.48,142.842105263158]
[86.925,4.777777777778,61.14,27.363157894737]
[3.275,-13.592694063927,61.62,-64.854736842105]
[40.875,0.315068493151,66.7,42.552631578947]
Press any key to continue . . . 

C doesn't do overloaded operators so I had to use the mangled name, and declare a large fixed array in the structure to reserve enough memory for fb dynamic arrays.

[Luxan]

A few more routines included, for data generation and Vector of Vectors print.

The whole structure of the matrices, and vectors, is determined by the
numerical values assigned to the array m_seq() .

Now using Matrix of matrices, and Vector of vectors, to do
Successive matrix vector multiplications with various dimensions for the
matrices and vectors.

Requires examination.
I have yet to compare these results with Maxima CAS .

Three pieces of code included.

------------------------( 1 )-----------------------

Code: Select all

' easy_matx.bas

' On my computer
' cd Downloads/fb/matrices_x/dodicat/merge_1/lib
'

' ' compile with: fbc -lib easy_matx.bas

'   Matrix math
'
'  1. dim as Matrix a
'  2. redim a.m( x1 to x2 , y1 to y2 )
'
'  3. dim as Matrix z(z1 to z2)  ?
'  4. redim z.m(w1 to w2)  ?
'
'     Matrix(x , y) ~= redim (0 to x-1, 0 to y-1) ?


' These declarations from :
'  http://www.rosettacode.org/wiki/Matrix_multiplication#FreeBASIC
' That content is available under
'  Content is available under GNU Free Documentation License 1.2 unless otherwise noted.
'
'   Using Geany, Build, Set Build Commands
'
'   Compile                fbc -w all "%f" 
'
'   Execute                "./%e" 
'
' Preliminary results checked using Maxima CAS
'
' also at, especially for vector by matrix calculations.
' https://keisan.casio.com/exec/system/15052033860538
'
' https://matrix.reshish.com/multCalculation.php
'
' https://elsenaju.eu/Calculator/matrix-vector-product.htm
'
'
type Matrix
    dim as double m( any , any )
    declare constructor ( )
    declare constructor ( byval x as uinteger , byval y as uinteger )
end type
 
constructor Matrix ( )
end constructor
 
constructor Matrix ( byval x as uinteger , byval y as uinteger )
    redim this.m( x - 1 , y - 1 )
end constructor


type Vector
dim as double v( any )
    declare constructor ( )
    declare constructor ( byval x as uinteger )
end type

constructor Vector ( )
end constructor

constructor Vector ( byval x as uinteger )
    redim this.v( x - 1 )
end constructor
'
' ====================================================================
'
Public operator * ( byref a as Matrix , byref b as Matrix ) as Matrix
'  Only applicable to square matrices or arrays.
' This routine from :
'  http://www.rosettacode.org/wiki/Matrix_multiplication#FreeBASIC
'
' Dodicat suggests only using: if ubound( a.m , 2 ) = ubound( b.m , 1 )
'
'
    dim as Matrix ret
    dim as uinteger i, j, k
    'if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
     if ubound( a.m , 2 ) = ubound( b.m , 1 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                for k = 0 to ubound( b.m , 1 )
                    ret.m( i , j ) += a.m( i , k ) * b.m( k , j )
                next k
            next j
        next i
    end if
    return ret
end operator
'
'
Public operator * ( byref a as Matrix , byref b as Vector ) as Vector
'
'   Vector multiplied by a matrix, maybe correct .
'
'
   dim as Vector ret
   dim as uinteger i, j, k
   if ubound( a.m , 2 ) = ubound( b.v , 1 )  then
        redim ret.v( ubound( a.m , 1 ) )
        for i = 0 to ubound( a.m , 1 )
            
                for k = 0 to ubound( b.v , 1 )
                    ret.v( i ) += a.m( i , k ) * b.v( k  )
                next k
            
        next i
    end if   
   return ret
end operator
'
Public operator - ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                ret.m( i , j ) = a.m( i , j ) - b.m( i , j )
            next j
        next i
    end if
    return ret
end operator
 '
Public operator + ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
               ret.m( i , j ) = a.m( i , j ) + b.m( i , j )
            next j
        next i
    end if
    return ret
end operator
'
Public operator / ( byref a as Matrix , byref b as Matrix ) as Matrix
    dim as Matrix ret
    dim as uinteger i, j, k
    dim as double x, y, z
    if ubound( a.m , 2 ) = ubound( b.m , 1 ) and ubound( a.m , 1 ) = ubound( b.m , 2 ) then
        redim ret.m( ubound( a.m , 1 ) , ubound( b.m , 2 ) )
        for i = 0 to ubound( a.m , 1 )
            for j = 0 to ubound( b.m , 2 )
                x = a.m( i , j )
                y = b.m( i , j )
                y = sgn(x)*10^32 ' assume y = 0
             if y <> 0 then z = x/y
                ret.m( i , j ) = z
            next j
        next i
    end if
    return ret
end operator
'
' ......................................................................
'
Public sub prt_z(z() as Matrix, idx as Integer)
'
'     Matrix of Matrix, index and print .
'
dim as integer i,j
'
'print
for i=0 to ubound(z(idx).m,1)
 for j=0 to ubound(z(idx).m,2) 
    print Using "###.####";z(idx).m(i,j);" ";
  next j
  print
next i
print
'  
end sub
'
' ......................................................................
'
Public sub prt_m(z as Matrix)
'
'  Print elements from rectangular matrix .
'
dim as integer i,j
'

print
print ubound(z.m,1);" ";ubound(z.m,2)
print
 for i=0 to ubound(z.m,1) 
  for j=0 to ubound(z.m,2)  
    print Using "###.####";z.m(i,j);" ";
  next j
  print
next i
print
'  
end sub
'
' ......................................................................
'
Public sub prt_v(z as Vector)
'
'  Print elements from vector .
'
dim as integer i
'
'print
for i=0 to ubound(z.v,1)
    
    print Using "###.####";z.v(i);" ";
  
  
next i
print
print
'  
end sub
'
' ......................................................................
'
Public sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
'     Multiply vector by matrix .
'
dim as integer i,j,ub1, ub2
dim as single x,y,s
    
    ub1=ubound(M1.m,1)
    ub2=ubound(M1.m,2)
'    
    redim V3.v(0 to ub1)
'
        for j=0 to ub1
             s=0
         for i=0 to ub2
             x=V1.v(i)
             y=M1.m(j,i)
             s=s+x*y
         next i
          V3.v(j)=s
      next j
 '
end sub
'
' ......................................................................
'
Public sub gen_z(z() as Matrix, idx as Integer)
'
'     Matrix of Matrix, index and generate values .
'
dim as integer i,j
'
Randomize
  for j=0 to ubound(z(idx).m,2) 
   for i=0 to ubound(z(idx).m,1)
      z(idx).m(i,j) = -1 + 2*rnd
  next i
next j

'  
end sub
'
' ......................................................................
'
Public sub gen_v(u() as Vector, idx as Integer)
'
'     Vector of Vector, index and generate values .
'
dim as integer i,j
Randomize
'
for i=0 to ubound(u(idx).v,1)
   
      u(idx).v(i) = -1 + 2*rnd
  
next i

'  
end sub
'
' ......................................................................
'
Public sub prtv_v(u() as Vector, idx as Integer)
'
'     Vector of Vector, index and print values .
'
dim as integer i,j
'
for i=0 to ubound(u(idx).v,1)
   print Using "###.####";u(idx).v(i)
next i
print 
'  
end sub
'
' ......................................................................
'

----------------------------(( 2 ))-------------------------

Code: Select all

' easy_matx.bi


#inclib "easy_matx"

'  For some reason, may need to duplicate types, but not constructors,
' in this *.bi

type Matrix
    dim as double m( any , any )
    declare constructor ( )
    declare constructor ( byval x as uinteger , byval y as uinteger )
end type
'
type Vector
dim as double v( any )
    declare constructor ( )
    declare constructor ( byval x as uinteger )
end type

'
'
' ---------------------------- operators -------------------------------
'
declare operator * ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator * ( byref a as Matrix , byref b as Vector ) as Vector

declare operator - ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator + ( byref a as Matrix , byref b as Matrix ) as Matrix
declare operator / ( byref a as Matrix , byref b as Matrix ) as Matrix

'   The rest from luxan, sciwiseg@gmail.com
'
'  Probably need to make this content
'  available under GNU Free Documentation License 1.2 
'
'
' ............................. procedures .............................
'
declare sub gen_z(z() as Matrix, idx as Integer)
declare sub gen_v(u() as Vector, idx as Integer)

declare sub prtv_v(u() as Vector, idx as Integer)

declare sub prt_z(z() as Matrix, idx as Integer)
declare sub prt_m(z as Matrix)
declare sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)

declare sub prt_v(z as Vector)
'

-----------------------------(((3)))-----------------------------------

Code: Select all

'
'  matx_smv.bas
'
'   Successive matrix vector multiplications .
'

'' compile with: fbc matx_test.bas


#include once "easy_matx.bi"

'
' ----------------------------------------------------------------------
'
dim as integer m_seq(0 to 4)={8,6,3,5}
dim as integer lms,nx,ny,i
lms=ubound(m_seq,1)
'print(lms)
'
'    Matrix of matrices
'
dim as Matrix z(1 to lms)
'
'   Vector of vectors
'
dim as Vector u(0 to lms)
'
'   Dimension matrices and vectors .
'
i=0
nx=m_seq(i)
u(i) = Vector(nx)
u(i).v(0)=1:u(i).v(1)=2:u(i).v(2)=3
u(i).v(3)=4:u(i).v(4)=5:u(i).v(5)=6
u(i).v(6)=7:u(i).v(7)=8
'
'prt_v(u(0))
'
for i=0 to lms-1
   nx=m_seq(i)
   ny=m_seq(i+1)
   z(i+1) = Matrix(ny , nx)
   u(i) = Vector(nx)
next i
'
'  ............... Assign values to input vector & matrices ............
'
gen_v(u(), 0)
print
print " Input vector, a matrix in another version ? "
print
prtv_v(u(), 0)
'
print   " matrix of matrices values "
print
'
for i=1 to lms
    gen_z(z() , i)
    prt_z(z() , i)
next i
print
'
print " sucessive vectors , matrices & vector x matrix calculations"
print " u(i) = z(i) * u(i - 1) "
print
'
print " u(0) "
prtv_v(u(), 0)
for i=1 to lms - 1
   print " z(";i;") "
   prt_z(z() , i)
   u(i) = z(i)*u(i-1)
   print " u(";i;") "
   prtv_v(u(), i)   
next i


'  end command acts as global destructor 
' of allocated memory ?

end
'
' ======================================================================
'

[fxm]

Run-time error, because at line 16 of the third file:
dim as integer m_seq(0 to 4)={8,6,3,5}
'm_seq(4)', the fifth element of 'm_seq()' is not initialized, therefore equal to '0'.

That induces a run-time error in the loop:
for i=0 to lms-1
nx=m_seq(i)
ny=m_seq(i+1)
z(i+1) = Matrix(ny , nx)
u(i) = Vector(nx)
next i
For 'i = lms-1'
=> 'z(3) = Matrix(0,5)'

I think that the correct code should be for example:
dim as integer m_seq(0 to 3)={8,6,3,5}
or:
dim as integer m_seq(0 to 4)={8,6,3,5,?}

[dodicat]

What is the sense of creating libraries to test/debug your WHOLE code at this stage?
You are adding another layer of complications.
My dll was only to test exporting one function (initmatrix) to C and pascal, which involved varargs, and I wasn't sure about it.
As it turns out it works well in these two languages.
As for any of the other languages listed in rosetta code, I haven't a clue whether it will work or not.

Libraries and .bi files -- suggest at the end of your project
fxm has the patience of a saint, but maybe he did the same as me, run easy_matx.bas plus your test.
Already you are getting tied up using a zero element for a matrix work.

I suggest use this option at the start of your code
#cmdline "-exx"
And use an ide option which doesn't close the console on a crash.
This will greatly help you along to your final library.
Personally I think I will make the operators from functions.
e.g.
operator *(M1 as matrix,m2 as matrix) as matrix
return matmult(M1,M2)
end function
Where all the work is done in matmult.
fb operators in dll's do not pass well to C or Pascal, although pascal can still do A*B.
C has only about 32 keywords, so these manoeuvres have to be different in C.
Sorry for being so irritable.