Hi dodicat
It is going take a day or two before I work the code on this topic into what I want to do,
which is:
1. Generate [two dimensional] matrices of arbitrary size; indexed through a dedicated array.
These matrices to be available for element wise operations and to retain values until
conclusion of other procedures.
As an example , generate matrices M[nx,ny] with index, through array, [0,1,2,3..] .
index1 <--- aidx[0]
M[8,6]
index2 <--- aidx[1]
M[6,4]
index3 <--- aidx[2]
M[4,7]
Equivalent to an indexing like M[idx,jdx](kdx)
where kdx is the index to the array that holds the matrices.
2. Perform various operations upon M[idx,jdx](kdx) , like addition, subtraction, multiplication;
other.
Similarly an indexed array of vectors is required.
Easy matrix
Re: Easy matrix
As I suggested before, why not an array of matrices, which index themselves by their array position.
(I use array of matrix 1 to whatever, but 0 to whatever would be OK)
(I use array of matrix 1 to whatever, but 0 to whatever would be OK)
Code: Select all
screen 20
Type matrix
As Double e(Any,Any)
End Type
Sub setsize(m As matrix,d1 As Long,d2 As Long)
Redim m.e(1 To d1,1 To d2)
End Sub
Sub show(z As matrix)
For n As Long=1 To Ubound(z.e,1)
For m As Long=1 To Ubound(z.e,2)
Print z.e(n,m);" ";
Next m
Print
Next n
End Sub
Dim As matrix z(1 To 5) 'ARRAY OF MATRICES
randomize
do
For n As Long=Lbound(z) To Ubound(z)
setsize(z(n),2+Rnd*4,2+Rnd*4)
Next
For n As Long =Lbound(z) To Ubound(z)
Print "matrix ";n
show(z(n))
Print
Next
For n As Long=Lbound(z) To Ubound(z)-1
For m As Long=n+1 To Ubound(z)
If Ubound(z(n).e,2)<>Ubound(z(m).e,1) Then
color 3
Print "matrix ";n;" and ";"matrix ";m; " won't multiply"
Else
color 15
Print "matrix ";n;" and ";"matrix ";m; " will multiply"
End If
Next
Next
color 15
print "Press a key or escape"
sleep
cls
loop until inkey=chr(27)
Re: Easy matrix
Example from: (How a function return can alone provide an array) thread
Code: Select all
#define lim 100
dim as double m(lim,lim)
m(0,1)=6 'rows
m(1,0)=8 'columns
'm(0,0)is spare
for r as long=1 to 6
for c as long=1 to 8
m(r,c)=int(rnd*50)
next c
next r
dim as double m2(lim,lim)
m2(0,1)=8 'rows
m2(1,0)=3 'columns
for r as long=1 to 8
for c as long=1 to 3
m2(r,c)=int(rnd*50)
next c
next r
Sub matmult(m1() As Double,m2() As Double,ans() As Double)
Dim rows As Long=m2(0,1)
Dim columns As Long=m1(1,0)
If columns<>rows Then Print "Can't do matmult":Return
Redim ans(lim,lim)
Dim rxc As Double
For r As Long=1 To rows
For c As Long=1 To columns
rxc=0
For k As Long = 1 To Ubound(m1,2)
rxc=rxc+m1(r,k)*m2(k,c)
Next k
ans(r,c)=rxc
Next c
Next r
ans(0,1)=m1(0,1)
ans(1,0)=m2(1,0)
End Sub
Sub show(A() As Double)
For rows As Long=1 To A(0,1)
For columns As Long=1 To A(1,0)
Print A(rows,columns);" ";
Next
Print
Next
End Sub
dim as double ans()
matmult(m(),m2(),ans())
show(m())
print "Multiply by"
show(m2())
print "Gives"
show(ans())
sleep
Re: Easy matrix
Using your code, and a little of my own, I'm able to observe what you've done; and how
it is relevant to what I want to do.
What you've done is fairly ingenious.
Merging the previous matrix procedures, and testing those with this, is the next step.
it is relevant to what I want to do.
What you've done is fairly ingenious.
Merging the previous matrix procedures, and testing those with this, is the next step.
Code: Select all
screen 20
Type matrix
As Double e(Any,Any)
End Type
Sub setsize(m As matrix,d1 As Long,d2 As Long)
Redim m.e(1 To d1,1 To d2)
End Sub
Sub show(z As matrix)
For n As Long=1 To Ubound(z.e,1)
For m As Long=1 To Ubound(z.e,2)
Print z.e(n,m);" ";
Next m
Print
Next n
End Sub
Dim As matrix z(1 To 5) 'ARRAY OF MATRICES
Dim as Integer n,nx,ny
n=1
nx=8
ny=6
setsize(z(n),nx,ny)
n=2
nx=6
ny=3
setsize(z(n),nx,ny)
'
n=1
Print "matrix ";n
show(z(n))
Print
'
n=2
Print "matrix ";n
show(z(n))
Print
'
sleep()
end
Re: Easy matrix
And what about an array of vectors.
Might arrays of matrices be of type Complex.
Might arrays of matrices be of type Complex.
Re: Easy matrix
Something I put together in the last few minutes, for array of vectors.
Code: Select all
' As I [ dodicat ] suggested before, why not an array of matrices,
' which index themselves by 'their array position.
' (I use array of matrix 1 to whatever, but 0 to whatever
' would be OK)
' modified for array of vectors .
screen 20
'
'
Type vector
As Double e(Any)
End Type
declare sub initial_V(v As vector)
declare Sub setsize_V(v As vector,d1 As Long)
declare Sub show_V(v As vector)
declare sub prt2v(v as vector, w as vector)
'
' ----------------------------------------------------------------------
'
Dim As vector v(1 To 5) 'ARRAY OF VECTORS
Dim as Integer n,nx,i ,j
n=1
nx=8
setsize_V(v(n),nx)
n=2
nx=6
setsize_V(v(n),nx)
'
n=1
initial_V(v(n))
Print "vector ";n
show_V(v(n))
Print
'
n=2
initial_V(v(n))
Print "vector ";n
show_V(v(n))
Print
'
'sleep()
n=2
i=3
j=1
print "v(";n;")(";i;")"
print ( v(n).e(i) )
n=1
i=5
j=3
print "v(";n;")(";i;")"
print ( v(n).e(i) )
print
prt2v(v(1), v(2))
'
sleep()
end
'
' ----------------------------------------------------------------------
'
sub prt2v(v as vector, w as vector)
'
' print element from different vectors
'
print "v(2)"
print ( v.e(2) )
print "w(3)"
print ( w.e(3) )
end sub
'
'
Sub setsize_V(v As vector,d1 As Long)
Redim v.e(1 To d1)
End Sub
'
Sub show_V(v As vector)
'
'
'
For n As Long=1 To Ubound(v.e,1)
Print v.e(n)
Next n
Print
End Sub
'
sub initial_V(v As vector)
dim as integer ix,lx
'
' sequential data for v( , )
'
lx=Ubound(v.e,1)
for ix=1 to lx
v.e(ix)=ix
next ix
'
end sub
'
'
Re: Easy matrix
I've been a little busy, and now have some code for you to examine and use.
Note that I decided to go with a Matrix of matrices,, for differently dimension-ed matrices indexing, as this uses the existing constructors.
There's a version without constructors or operators, that I might eventually complete.
Note that I decided to go with a Matrix of matrices,, for differently dimension-ed matrices indexing, as this uses the existing constructors.
There's a version without constructors or operators, that I might eventually complete.
Code: Select all
' 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
'
'
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
'
' ---------------------------- operators -------------------------------
'
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
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)
'
' ----------------------------------------------------------------------
'
' 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)
'
'
prt_z(z() , 2)
prt_z(z() , 5)
end
'
' ====================================================================
'
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
'
'
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
'
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
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
' Multiply vector by matrix .
'
dim as integer i,j,k,ub1, ub2,vb1
dim as single x,y,s
'vb1=ubound(M1.m,1)
ub1=ubound(M1.m,1)
ub2=ubound(M1.m,2)
'
redim V3.v(0 to ub2+1)
'
for j=0 to ub2+1
s=0
for i=0 to ub1
x=V1.v(i)
y=M1.m(j,i)
s=s+x*y
next i
V3.v(j)=s
next j
'
end sub
'
' ......................................................................
'
Re: Easy matrix
1)
2)
Instead of:
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
You can define another overload '*' operator:
operator * ( byref a as Matrix , byref b as Vector ) as Vector
Code: Select all
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
' Multiply vector by matrix .
'
dim as integer i,j,k,ub1, ub2,vb1
dim as single x,y,s
'vb1=ubound(M1.m,1)
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)
Instead of:
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
You can define another overload '*' operator:
operator * ( byref a as Matrix , byref b as Vector ) as Vector
Re: Easy matrix
Apart from the vector multiplication, is there anything else I should be aware of .
From Maxima CAS, this is what I get when I multiply a vector by a matrix.
M1,V1
matrix(
[14],
[32],
[50],
[18],
[13]
)
This is also what I got when I ran the FreeBasic code.
Perhaps, I can treat a vector as a single row, or column, matrix; might require redim and
extending with zero elements.
From Maxima CAS, this is what I get when I multiply a vector by a matrix.
Code: Select all
kill(all)$
ratprint:false$
print(" M1 ")$
M1:matrix([1,2,3],[4,5,6],[7,8,9],[7,1,3],[5,1,2]);
print(" ")$
V1:matrix([1,2,3]);
print(" M1.V1 ")$
M1.V1;
matrix(
[14],
[32],
[50],
[18],
[13]
)
This is also what I got when I ran the FreeBasic code.
Perhaps, I can treat a vector as a single row, or column, matrix; might require redim and
extending with zero elements.
Re: Easy matrix
Using vector with matrix is only a convenience for me, which is OK for solving linear equations and using in inverse, regression e.t.c.
But bare matrix work is more logical using matrices only.
But bare matrix work is more logical using matrices only.
Code: Select all
Dim Shared As boolean ExtremePivot
const decplaces=6
ExtremePivot=true
#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))
Type vector
Dim As Double element(Any)
Declare Operator Cast() As String
End Type
Type matrix
Dim As Double element(Any,Any)
Declare Operator Cast() As String
Declare Function inverse() As matrix
Declare Function GaussJordan(As vector) As vector
Declare Function determinant() As Double
Declare Function transpose() As matrix
Declare Function unit(r As Long) As matrix
Declare Constructor(As Long,As Long)
declare constructor
Declare Function vandermonde(() As Double) As matrix
declare sub set(as long,as long,as double)
End Type
sub matrix.set(r as long,c as long,x as double)
this.element(r,c)=x
end sub
constructor matrix
end constructor
Constructor matrix (r As Long,c As Long)
Redim element(1 To r,1 To c)
End Constructor
Function matrix.unit(r As Long) As matrix
Dim As matrix m=Type(r,r)
For z As Long=1 To r
m.element(z,z)=1
Next
Return m
End Function
Operator vector.cast() As String 'for printing
Dim As String ans
ans= "Solution:"+Chr(10)
For n As Integer=1 To Ubound(this.element)
ans=ans+ Str(Cround(element(n),decplaces))+Chr(10)
Next n
Operator =ans
End Operator
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
Function matrix.transpose() As matrix
Dim As matrix b
Redim b.element(1 To Ubound(this.element,2),1 To Ubound(this.element,1))
For i As Long=1 To Ubound(this.element,1)
For j As Long=1 To Ubound(this.element,2)
b.element(j,i)=this.element(i,j)
Next
Next
Return b
End Function
Operator *(m1 As matrix,m2 As matrix) As matrix
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
Operator *(m1 As matrix,m2 As vector) As vector
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element) Then
Print "Can't do multiply"
Exit Operator
End If
Dim As vector ans
Redim ans.element(rows)
Dim rxc As Double
For r As Integer=1 To rows
rxc=0
For k As Integer = 1 To Ubound(m1.element,2)
rxc=rxc+m1.element(r,k)*m2.element(k)
Next k
ans.element(r)=rxc
Next r
Operator= ans
End Operator
Operator +(m1 As matrix,m2 As matrix) As matrix
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,1)<>Ubound(m2.element,1) or Ubound(m1.element,2)<>Ubound(m2.element,2) Then
Print "Can't do plus"
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
ans.element(r,c)=m1.element(r,c)+m2.element(r,c)
Next c
Next r
Operator= ans
End Operator
Operator -(m1 As matrix,m2 As matrix) As matrix
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,1)<>Ubound(m2.element,1) or Ubound(m1.element,2)<>Ubound(m2.element,2) Then
Print "Can't do minus"
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
ans.element(r,c)=m1.element(r,c)-m2.element(r,c)
Next c
Next r
Operator= ans
End Operator
Function matrix.GaussJordan(rhs As vector) As vector
Dim As Integer n=Ubound(rhs.element)
Dim As vector ans=rhs,r=rhs
Dim As matrix b=This
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b.element(p1,num))<Abs(b.element(p2,num)) Then
Swap r.element(p1),r.element(p2)
For g As Integer=1 To n
Swap b.element(p1,g),b.element(p2,g)
Next g
End If
Next p2
Next p1
#endmacro
If ExtremePivot=false Then :pivot(1): End If
For k As Integer=1 To n-1
If ExtremePivot=true Then: pivot(k) :End If
For row As Integer =k To n-1
If b.element(row+1,k)=0 Then Exit For
Var f=b.element(k,k)/b.element(row+1,k)
r.element(row+1)=r.element(row+1)*f-r.element(k)
For g As Integer=1 To n
b.element((row+1),g)=b.element((row+1),g)*f-b.element(k,g)
Next g
Next row
Next k
'back substitute
For z As Integer=n To 1 Step -1
ans.element(z)=r.element(z)/b.element(z,z)
For j As Integer = n To z+1 Step -1
ans.element(z)=ans.element(z)-(b.element(z,j)*ans.element(j)/b.element(z,z))
Next j
Next z
Function = ans
End Function
Function matrix.inverse() As matrix
If Ubound(element,1)<>Ubound(element,2) Then
Print "Can't do inverse"
return unit(0)
end if
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2)
Dim As matrix ans
Dim As vector temp,null_
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1)
Redim ans.element(1 To ub1,1 To ub2)
For a As Integer=1 To ub1
temp=null_
temp.element(a)=1
temp=GaussJordan(temp)
For b As Integer=1 To ub1
ans.element(b,a)=temp.element(b)
Next b
Next a
Return ans
End Function
Function matrix.determinant() As Double
If Ubound(element,1)<>Ubound(element,2) Then
Print "Can't do determinant"
Exit function
end if
Dim As Integer n=Ubound(this.element,1),sign=1
Dim As Double det=1,s=1
Dim As Double b(1 To n,1 To n)
For c As Integer=1 To n
For d As Integer=1 To n
b(c,d)=this.element(c,d)
Next d
Next c
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b(p1,num))<Abs(b(p2,num)) Then
sign=-sign
For g As Integer=1 To n
Swap b(p1,g),b(p2,g)
Next g
End If
Next p2
Next p1
#endmacro
If ExtremePivot=false Then :pivot(1): End If
For k As Integer=1 To n-1
If ExtremePivot =true Then: pivot(k) :End If
For row As Integer =k To n-1
If b(row+1,k)=0 Then Exit For
Var f=b(k,k)/b(row+1,k)
s=s*f
For g As Integer=1 To n
'divide by f here seems ok
b((row+1),g)=(b((row+1),g)*f-b(k,g))/f
Next g
Next row
Next k
For z As Integer=1 To n
det=det*b(z,z)
Next z
Return sign*det
End Function
Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
Var n=Ubound(x_values)
Dim As matrix Mat
Redim mat.element(1 To n,1 To n)
Dim As vector rhs,ans
Redim rhs.element(1 To n)
Redim ans.element(1 To n)
Redim p(0):Redim p(1 To n)
For a As Integer=1 To n
rhs.element(a)=y_values(a)
For b As Integer=1 To n
mat.element(a,b)=x_values(a)^(b-1)
Next b
Next a
'Solve the linear equations
ans=mat.GaussJordan(rhs)
For z As Integer=1 To n:p(z)=ans.element(z):Next
End Sub
Sub setmatrix(m As matrix,d() As Double)
m=Type(Ubound(d,1),Ubound(d,2))
For n1 As Long=1 To Ubound(d,1)
For n2 As Long=1 To Ubound(d,2)
m.element(n1,n2)=d(n1,n2)
Next
Next
End Sub
'vandermode of x
Function vandermonde(x_values() As Double,w As Long) As matrix
Dim As matrix mat
Var n=Ubound(x_values)
Redim mat.element(1 To n,1 To w)
For a As Integer=1 To n
For b As Integer=1 To w
mat.element(a,b)=x_values(a)^(b-1)
Next b
Next a
Return mat
End Function
Sub regress(x_values() As Double,y_values() As Double,ans() As Double,n As Long)
Redim ans(1 To Ubound(x_values))
Dim As matrix m1= vandermonde(x_values(),n)
Dim As matrix T=m1.transpose
Dim As vector y
Redim y.element(1 To Ubound(ans))
For n As Long=1 To Ubound(y_values)
y.element(n)=y_values(n)
Next n
Dim As vector result=(((T*m1).inverse)*T)*y
Redim Preserve ans(1 To n)
For n As Long=1 To Ubound(ans)
ans(n)=result.element(n)
Next n
End Sub
'Evaluate a polynomial at x
Function polyeval(Coefficients() As Double,Byval x As Double) As Double
Dim As Double acc
For i As Long=Ubound(Coefficients) To Lbound(Coefficients) Step -1
acc=acc*x+Coefficients(i)
Next i
Return acc
End Function
'=========== SET UP A SYSTEM ======================
'=========== SET UP A SYSTEM ======================
Dim As Double d1(1 To 5,1 To 3)={ {1,2,3},_
{4,5,6}, _
{7,8,9}, _
{7,1,3}, _
{5,1,2} }
Dim As matrix a1,a2
setmatrix(a1,d1())
print "matrix 1"
Print a1
Dim As Double d2(1 To 3,1 To 1)={ {1}, _
{2}, _
{3}}
setmatrix(a2,d2())
print "matrix 2"
print a2
print "matrix 1 * matrix 2"
print a1*a2
print "----------------------------"
Dim As Double d3(1 To 1,1 To 3)={ {1,2,3} }
setmatrix(a1,d3())
print "matrix 1"
print a2
print "matrix 2"
print a1
print "matrix 1 * matrix 2"
print a2*a1
print "matrix 2 * matrix 1"
print a1*a2
sleep
Re: Easy matrix
Vector by matrix multiplication, using operator overloading; again.
I wonder whether operator overloading means a slower program.
I now feel a need to examine vector by matrix multiplication more thou roughly, perhaps
there are interactive website that I might use.
I wonder whether operator overloading means a slower program.
I now feel a need to examine vector by matrix multiplication more thou roughly, perhaps
there are interactive website that I might use.
Code: Select all
' 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
'
'
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
'
' ---------------------------- 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)
'
' ----------------------------------------------------------------------
'
' 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
'
'
z(1)=a
z(2)=M1
z(3)=M2
'z(4)=M2
'
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)
'
'
prt_z(z() , 2)
prt_z(z() , 5)
end
'
' ====================================================================
'
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
'
'
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
'
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
'
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
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
' Multiply vector by matrix .
'
dim as integer i,j,k,ub1, ub2,vb1
dim as single x,y,s
'vb1=ubound(M1.m,1)
ub1=ubound(M1.m,1)
ub2=ubound(M1.m,2)
'
redim V3.v(0 to ub2+1)
'
for j=0 to ub2+1
s=0
for i=0 to ub1
x=V1.v(i)
y=M1.m(j,i)
s=s+x*y
next i
V3.v(j)=s
next j
'
end sub
'
' ......................................................................
'
Re: Easy matrix
Vector by matrix calculations validated at these websites.
https://keisan.casio.com/exec/system/15052033860538
https://matrix.reshish.com/multCalculation.php
https://elsenaju.eu/Calculator/matrix-v ... roduct.htm
https://keisan.casio.com/exec/system/15052033860538
https://matrix.reshish.com/multCalculation.php
https://elsenaju.eu/Calculator/matrix-v ... roduct.htm
Re: Easy matrix
These are linear equation type calculations, the vector type (and vector multiply) is handy for these.
Like everything else, get something to do a task easily.
But for mathematical representation, only matrices are better.
They are a little more difficult to set up with that extra dimension 1 to 1 for a vector type, but they shouldn't let you down because they follow the rules of matrices.
Like everything else, get something to do a task easily.
But for mathematical representation, only matrices are better.
They are a little more difficult to set up with that extra dimension 1 to 1 for a vector type, but they shouldn't let you down because they follow the rules of matrices.
Re: Easy matrix
The execution time is not penalized because the overload resolution is done at the compile time, only depending on the declared type of the arguments passed to the overload procedure.
In your last code, you did not correct the 'Vect_x_Matrix()' procedure that always induces a runtime error (seen with the '-exx' compilation option).
See my proposed correction in my previous posts.
Re: Easy matrix
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.
Vector as a matrix is the next to do.
As requested, here's there updated code, utilizing your suggestions.
Incorporated into the rest of the code like this, is it alright with you.
Vector as a matrix is the next to do.
Code: Select all
' 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
'
' ---------------------------- 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)
'
' ----------------------------------------------------------------------
'
' 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)
'
'
prt_z(z() , 2)
prt_z(z() , 5)
end
'
' ====================================================================
'
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
'
'
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
'
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
'
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
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
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
'
' ......................................................................
'
sub Vect_x_Matrix(M1 as Matrix, V1 as Vector, V3 as Vector)
'
' Multiply vector by matrix .
'
dim as integer i,j,k,ub1, ub2,vb1
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