After installing fbmath, with all of its fancy matrix routines,
I encounter a difficulty in finding therein a routine for matrix multiplication, addition,subtraction,division,
and vector by matrix multiplication .
Attempted the installation of the fb extended math library, became complicated
, and ineffective , rather quickly.
This illustrates the need for better and simpler code for what are really
basic routines, that are only dependent upon the standard libraries that are
installed with FreeBasic.
Easy matrix
Re: Easy matrix
Normally I prefer simple arrays of double for matrix routines.
However I made up a few routines years ago, under a structure (type).
EDIT
Put some dimension guards in.
However I made up a few routines years ago, under a structure (type).
EDIT
Put some dimension guards in.
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=1,As Long=1)
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 (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=Ubound(this.element,2) Then comma="" Else comma=","
Var q=Str(Cround(element(a,b),decplaces))
ans=ans+q+comma'+gap(q)
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"
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 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 ======================
Dim As Double d1(1 To 5,1 To 5)={ {2,7,9,0,1},_
{0,7,-6,2,7}, _
{9,9,1,-5,2}, _
{0,9,0,9,8}, _
{6,7,7,8,9} }
Dim As matrix a1,a2
setmatrix(a1,d1())
Print a1
Dim As Double d2(1 To 5,1 To 5)={ {2,7,9,0,1},_
{0,7,-6,2,7}, _
{9,9,1,-5,2}, _
{0,9,0,9,8}, _
{6,7,7,8,9} }
setmatrix(a2,d2())
Print a2
Print "PLUS (the two at top)"
Print a1+a2
Print "MINUS (the two at top)"
Print a1-a2
Print "inverse of PLUS the two at top"
Print(a1+a2).inverse
Print "Multiply top two"
Print a1*a2
Print "Transpose of Multiply top two"
Print (a1*a2).transpose
print "determinant of first matrix at the top"
print a1.determinant
Print "Press a key to continue . . ."
Sleep
Dim As Double x(1 To ...)={50,70,90,200,350,500,600}
Dim As Double y(1 To ...)={150,200,300,200,90,20,400}
Redim As Double ans()
regress(x(),y(),ans(),4) '4 terms = (power 3)
Screen 20
For n As Long=1 To Ubound(x)
Circle(x(n),y(n)),3,,,,,f
Next
For x As Long=0 To 1000 Step 1
Var f=polyeval(ans(),x)
Circle(x,f),1,2
Next
Locate 25
Color 2
Print "Regression curve through white dots"
For n As Long=Lbound(ans) To Ubound(ans)
If n=Lbound(ans) Then Print "Constant",ans(n) _
Else Print "x^";n-1,ans(n)
Next
Sleep
Last edited by dodicat on Jul 26, 2022 23:10, edited 1 time in total.
Re: Easy matrix
I've been to the Rosetta Code site and found some Freebasic code there, for matrix multiplication.
Need to check copyright conditions for that.
I added other routines for matrix addition, subtraction and division; based upon the sample code there.
I don't mind sharing my code, under a LGPL, with the rest of the community; just busy at the moment.
Need to check copyright conditions for that.
I added other routines for matrix addition, subtraction and division; based upon the sample code there.
I don't mind sharing my code, under a LGPL, with the rest of the community; just busy at the moment.
Re: Easy matrix
Technically there is no such thing as matrix division.
I left it as an exercise to concoct a multiplication by an inverse in about three lines, and call it a division operator.
The only copyright I would suggest for my code is just copy it right.
I left it as an exercise to concoct a multiplication by an inverse in about three lines, and call it a division operator.
The only copyright I would suggest for my code is just copy it right.
Re: Easy matrix
The division I did was just element by corresponding element in the other matrix.
You mean just copy it without any type of copyright.
Now I want to construct two dimensional arrays of various sizes, the name of each array to
be indexed by a more global variable; how might I do this.
You mean just copy it without any type of copyright.
Now I want to construct two dimensional arrays of various sizes, the name of each array to
be indexed by a more global variable; how might I do this.
Re: Easy matrix
I have put some guards in the code to check the dimensions.
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=1,As Long=1)
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 (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=Ubound(this.element,2) Then comma="" Else comma=","
Var q=Str(Cround(element(a,b),decplaces))
ans=ans+q+comma'+gap(q)
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"
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 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 ======================
Dim As Double d1(1 To 3,1 To 5)={ {2,7,9,0,1},_
{0,7,-6,2,7}, _
{9,9,1,-5,2}}
Dim As matrix a1,a2
setmatrix(a1,d1())
print !"a1\n"
Print a1
Dim As Double d2(1 To 5,1 To 3)={ {2,7,9},_
{0,7,-6}, _
{9,9,1}, _
{0,9,0}, _
{6,7,7} }
setmatrix(a2,d2())
print !"a2\n"
print a2
print !"a1*a2\n"
print a1*a2
print !"a2*a1\n"
print a2*a1
dim as matrix m=matrix(2,3)
m.set(1,1,rnd)
m.set(1,2,rnd)
m.set(1,3,rnd)
m.set(2,1,rnd)
m.set(2,2,rnd)
m.set(2,3,rnd)
print !"m\n"
print m
print !"m.transpose\n"
print m.transpose
print !"(m.transpose)*m\n"
print (m.transpose)*m
print !"m*(m.transpose)\n"
print m*(m.transpose)
print "Checking guards"
print m.determinant
print m.inverse
print m+a1
sleep
Re: Easy matrix
I want to dimension an arbitrary number of arrays, each with independent sizes for
their dimensions.
Then I need to index each array, as a separate entity.
I'm able to get a pointer to an array that I re dimension within a function.
To get the correct pointers to all of the arrays, adding in an offset of (n+m)*sizeof(var)
is suggested.
These pointers to be stored in a separate array.
their dimensions.
Then I need to index each array, as a separate entity.
I'm able to get a pointer to an array that I re dimension within a function.
To get the correct pointers to all of the arrays, adding in an offset of (n+m)*sizeof(var)
is suggested.
These pointers to be stored in a separate array.
Re: Easy matrix
Do your arrays of different sizes nevertheless all have the same dimension (same number of indexes) ?
Re: Easy matrix
I'm using 2 dimensional arrays and 1 dimensional arrays, these are destined to be used
as matrices and vectors.
Also looking at what others have done with constructors, something I'm not familiar with.
as matrices and vectors.
Also looking at what others have done with constructors, something I'm not familiar with.
Re: Easy matrix
I propose a structure for indexed vectors and another for indexed matrices:
Code: Select all
Type Vector
Dim As Integer element(Any)
Declare Sub Resize(Byval n As Integer)
Declare Property Ubound() As Integer
End Type
Sub Vector.Resize(Byval n As Integer)
Redim Preserve This.element(n - 1)
End Sub
Property Vector.Ubound() As Integer
Return ..Ubound(This.element)
End property
'----------------------------------------------
Dim As Vector indexedVector()
' Creating a first vector (index 0)
Redim Preserve indexedVector(Ubound(indexedVector) + 1)
' Sizing this first vector (index 0)
indexedVector(0).Resize(2)
' Filling this first vector (index 0)
For I As Integer = 0 To indexedVector(0).Ubound
indexedVector(0).element(I) = I + 1
Next I
' Adding a seconf vector (index 1)
Redim Preserve indexedVector(Ubound(indexedVector) + 1)
' Sizing this second vector (index 1)
indexedVector(1).Resize(3)
' Filling this second vector (index 1)
For I As Integer = 0 To indexedVector(1).Ubound
indexedVector(1).element(I) = I + 1
Next I
' Printing all indexed vectors
For K As Integer = 0 To Ubound(indexedVector)
For I As Integer = 0 To indexedVector(K).Ubound
Print indexedVector(K).element(I),
Next I
Print
Next K
Sleep
Code: Select all
Type Matrix
Dim As Integer element(Any, Any)
Declare Sub Resize(Byval n1 As Integer, Byval n2 As Integer)
Declare Property Ubound(Byval n As Integer) As Integer
End Type
Sub Matrix.Resize(Byval n1 As Integer, Byval n2 As Integer)
Redim Preserve This.element(n1 - 1, n2 - 1)
End Sub
Property Matrix.Ubound(Byval n As Integer) As Integer
Return ..Ubound(This.element, n)
End property
'---------------------------------------------------------------
Dim As Matrix indexedMatrix()
' Creating a first matrix (index 0)
Redim Preserve indexedMatrix(Ubound(indexedMatrix) + 1)
' Sizing this first matrix (index 0)
indexedMatrix(0).Resize(2, 3)
' Filling this first matrix (index 0)
For I As Integer = 0 To indexedMatrix(0).Ubound(1)
For J As Integer = 0 To indexedMatrix(0).Ubound(2)
indexedMatrix(0).element(I, J) = I * 10 + J + 1
Next J
Next I
' Adding a seconf matrix (index 1)
Redim Preserve indexedMatrix(Ubound(indexedMatrix) + 1)
' Sizing this second matrix (index 1)
indexedMatrix(1).Resize(4, 5)
' Filling this second matrix (index 1)
For I As Integer = 0 To indexedMatrix(1).Ubound(1)
For J As Integer = 0 To indexedMatrix(1).Ubound(2)
indexedMatrix(1).element(I, J) = I * 10 + J + 1
Next J
Next I
' Printing all indexed matrices
For K As Integer = 0 To Ubound(indexedMatrix)
For I As Integer = 0 To indexedMatrix(K).Ubound(1)
For J As Integer = 0 To indexedMatrix(K).Ubound(2)
Print indexedMatrix(K).element(I, J),
Next J
Print
Next I
Print
Next K
Sleep
Re: Easy matrix
Looks more than promising.
I'm imagining how I might incorporate this into my other Freebasic code,
seems straight forward.
I'm imagining how I might incorporate this into my other Freebasic code,
seems straight forward.
Re: Easy matrix
LUXAN
You can have a simple array of matrices.
(Why do you need pointers to a matrix?)
In my previous matrix posts to you, could you add the empty constructor, as in the above, and
Declare Constructor(As Long,As Long) as the other constructor.
Also in the inverse
Print "Can't do inverse"
return unit(0)
end if
(return 0 will give an error, as it should have done without the empty constructor)
You can have a simple array of matrices.
(Why do you need pointers to a matrix?)
Code: Select all
const decplaces=3
#define gap(x) string(5*decplaces-len(x)," ")
#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 matrix
Dim As Double element(Any,Any)
Declare Operator Cast() As String
Declare Constructor(As Long,As Long)
declare constructor
declare sub set(as long,as long,as double)
declare operator +=(as matrix)
End Type
sub matrix.set(r as long,c as long,x as double)
this.element(r,c)=x
end sub
Operator matrix.cast() As String 'for printing
Dim As String ans,comma
For a As Integer=1 To Ubound(this.element,1)
For b As Integer=1 To Ubound(this.element,2)
If b=Ubound(this.element,2) Then comma="" Else comma=","
Var q=Str(Cround(element(a,b),decplaces))
ans=ans+q+comma+gap(q)
Next b
ans=ans+Chr(10)
Next a
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)
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 matrix.+=(m as matrix)
this= this+m
end operator
Constructor matrix (r As Long,c As Long)
Redim element(1 To r,1 To c)
End Constructor
constructor matrix
end constructor
randomize
dim as long lim=1000
dim as matrix m(1 to lim) 'ARRAY OF MATRICES
for n as long=1 to lim
m(n)=matrix(2,2)
for z1 as long=1 to 2
for z2 as long=1 to 2
m(n).set(z1,z2,(200*(rnd-rnd)))
next z2
next z1
if n < 6 or n>lim-5 then
print "matrix ";n
print m(n)
else
locate csrlin-1
print "------------------"
end if
next n
print
dim as matrix tot=matrix(2,2)
for n as long=1 to lim
tot+=m(n)
next
print "All added up"
print tot
sleep
Declare Constructor(As Long,As Long) as the other constructor.
Also in the inverse
Print "Can't do inverse"
return unit(0)
end if
(return 0 will give an error, as it should have done without the empty constructor)
Re: Easy matrix
Hi dodicat
Your approach is interesting, I'm learning new concepts.
Does the use of a constructor and elaborate types mean a longer run time.
The arrays I intend to use might be of type string, most likely though they'll be of type
single or double; these are of differing dimensions, your code might accommodate that.
The array of matrices approach is what I did with an entirely different language.
Your approach is interesting, I'm learning new concepts.
Does the use of a constructor and elaborate types mean a longer run time.
The arrays I intend to use might be of type string, most likely though they'll be of type
single or double; these are of differing dimensions, your code might accommodate that.
The array of matrices approach is what I did with an entirely different language.
Re: Easy matrix
In this above application case yes, because for each indexed matrix, 2 objects are built, one of which is temporary, then this temporary object is destroyed:
- 'dim as matrix m(1 to lim)'
- each all objects are default constructed ('dim as matrix') with a non-sized member array,
- for each of such objects, a temporary object is built with a sized member array ('matrix(2,2)'), this temporary object is then copied into the initial object ('m(n) = ...'), then this temporary object is destroyed.
Re: Easy matrix
Example of using arrays only for matrices.
(linear equations using determinants only)
(linear equations using determinants only)
Code: Select all
Dim Shared As boolean ExtremePivot
const decplaces=20
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))
'don't use cround here
Function determinant(matrix() As Double) As Double
Dim As Long n=Ubound(matrix,1),sign=1
Dim As Double det=1,s=1
Dim As Double b(1 To n,1 To n)
For c As Long=1 To n
For d As Long=1 To n
b(c,d)=matrix(c,d)
Next d
Next c
#macro pivot(num)
For p1 As Long = num To n - 1
For p2 As Long = p1 + 1 To n
If Abs(b(p1,num))<Abs(b(p2,num)) Then
sign=-sign
For g As Long=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 Long=1 To n-1
If ExtremePivot=true Then: pivot(k) :End If
For row As Long =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 Long=1 To n
b((row+1),g)=(b((row+1),g)*f-b(k,g))/f
Next g
Next row
Next k
For z As Long=1 To n
det=det*b(z,z)
Next z
Return sign*det
End Function
'CRAMER COLUMN SWAPS
Sub swapcolumn(m() As Double,c() As Double,_new() As Double,column As Long)
Redim _new(1 To Ubound(m,1),1 To Ubound(m,1))
For x As Long=1 To Ubound(m,1)
For y As Long=1 To Ubound(m,1)
_new(x,y)=m(x,y)
Next y
Next x
For z As Long=1 To Ubound(m,1)
_new(z,column)=c(z,1)
Next z
End Sub
Sub solve(mat() As Double,rhs() As Double,_out() As Double)
If Ubound(mat,1) <> Ubound(rhs,1) Then Print "Can't do solve":Return
Redim _out(Lbound(mat,1) To Ubound(mat,1),1 To 1)
Redim As Double result(Lbound(mat,1) To Ubound(mat,1),Lbound(mat,1) To Ubound(mat,1))
Dim As Double maindeterminant=determinant(mat())
If Abs(maindeterminant) < 1e-12 Then Print "singular":Return
For column As Long=1 To Ubound(mat,1)
swapcolumn(mat(),rhs(),result(),column)
_out(column,1)= determinant(result())/maindeterminant
Next
End Sub
Sub matmult(m1() As Double,m2() As Double,ans() As Double)
Dim rows As Long=Ubound(m1,1)
Dim columns As Long=Ubound(m2,2)
If Ubound(m1,2)<>Ubound(m2,1) Then Print "Can't do matmult":Return
Redim ans(1 To rows,1 To columns)
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
End Sub
Sub show(A() As Double)
For n As Long=1 To Ubound(A,1)
For m As Long=1 To Ubound(A,2)
' Print Cround(A(n,m),decplaces);" ";
Print A(n,m);" ";
Next
Print
Next
End Sub
Dim As Double MainMat(1 To ...,1 To ...)={{2,-1,5,1}, _
{3,2,2,-6}, _
{1,3,3,-1}, _
{5,-2,-3,3}}
Dim As Double rhs(1 To ...,1 To 1)= {{-3}, _
{-32}, _
{-47}, _
{49}}
Print "Main matrix"
show(MainMat())
Print "Right hand side vector"
show(rhs())
Redim ans() As Double
solve(MainMat(),rhs(),ans())
Print "Solution Main Matrix * (something) = Right hand side vector"
show(ans())
Print
Redim As Double res()
matmult(MainMat(),ans(),res())
Print "Check the right hand side vector is returned"
print "i.e. Matrix * (solution) = Right hand side vector"
show(res())
Print
Dim As Double d1(1 To 3,1 To 5)={ {2,7,9,0,1},_
{0,7,-6,2,7}, _
{9,9,1,-5,2}}
Dim As Double d2(1 To 5,1 To 3)={ {2,7,9},_
{0,7,-6}, _
{9,9,1}, _
{0,9,0}, _
{6,7,7} }
Print "Matrix 1"
show(d1())
Print "Matrix 2"
show(d2())
print
Redim As Double product()
matmult(d1(),d2(),product())
Print "Product Matrix 1 x Matrix 2"
show(product())
print
Print
erase(product)
matmult(d2(),d1(),product())
Print "Product Matrix 2 x Matrix 1"
show(product())
Print
Dim As Double a(1 To ...,1 To 1)= {{-3}, _
{-32}, _
{-47.5}, _
{49}, _
{-5.5}, _
{-1}, _
{10}}
Print "matrix a"
show(a())
Print
Dim As Double b(1 To 1,1 To ...)= {{-2,-7.7,21.1,2,0,6,5}}
Print "matrix b"
show(b())
Print
Redim As Double c()
matmult(a(),b(),c())
Print "a*b"
show(c())
Print
erase(c)
matmult(b(),a(),c())
Print "b*a"
show(c())
Sleep