classes add getval pow svd bcopy ident printf symmeig c inverse resize transpose exp mulm scanf x fprint muls setcol zero getcol mulv setdiag getdiag ncol setrow getrow nrow solv
mobj = new Matrix(nrow, ncol)mobj = new Matrix(nrow, ncol, type)Individual element values are assigned and evaluated using the syntax:
which may appear anywhere in an expression or on the left hand side of an assignment statement. irow can range from 0 to m.nrow-1 and icol ranges from 0 to m.ncol-1 . ( See x )m.x[irow][icol]
When possible, Matrix methods returning a Matrix use the form, mobj = m.f(args, [mout]), where mobj is a newly constructed matrix (m is unchanged) unless the optional last mout argument is present, in which case the return value is mout and mout is used to store the matrix. This style seems most efficient to me since many matrix operations cannot be done in-situ. Exceptions to this rule, eg m.zero(), are noted in the individual methods.
Similarly, Matrix methods returning a Vector use the form, vobj = m.f(args, [vout]), where vobj is a newly constructed Vector unless the optional last vout is present for storage of the vector elements. Use of vout is extremely useful in those cases where the vector is graphed since the result of the matrix operation does not invalidate the pointers to the vector elements.
Note that the return value allows these operations to be used as members of a filter chain or arguments to other functions.
By default, a new Matrix is of type MFULL (= 1) and allocates storage for
all nrow*ncol elements. Scaffolding is in place for matrices of storage
type MSPARSE (=2) and MBAND (=3) but not many methods have been interfaced
to the meschach library at this time. If a method is called on a matrix type
whose method has not been implemented, an error message will be printed.
It is intended that implemented methods will be transparent to the user, eg
m*x=b (x = m.solv(b) ) will solve the linear system
regardless of the type of m and
v1 = m*v2 (v1 = m.mulv(v2) ) will perform the vector multiplication.
Matrix is implemented using the
meschach c library by David E. Stewart
(discovered at
http://www.netlib.org/c/index.html) which contains a large collection
of routines for sparse, banded, and full matrices. Many of the useful
routines have not
been interfaced with the hoc interpreter but can be easily added on request
or you can add it yourself
by analogy with the code in nrn/src/ivoc/(matrix.c ocmatrix.[ch])
At this time the MFULL matrix type is complete enough to do useful work
and MSPARSE can be used to multiply a matrix by a vector and solve
Mx=b.
Matrix
val = m.x[irow][icol]m.x[irow][icol] = valexpr(m.x[irow][icol]) &m.x[irow][icol]which may appear anywhere in an expression or on the left hand side of an assignment statement. irow can range from 0 to m.nrow-1 and icol ranges from 0 to m.ncol-1 .m.x[irow][icol]
For functions that require the address of a double value one may write
but one must be on guard for the case in which matrix storage is freed while another object holds a pointer to one of its elements. (Matrix does not currently notify the interpreter when storage has been freed.)&m.x[irow][icol]
For sparse matrices an invocation of x[i][j] will create it if the element does not exist. Therefore if you wish to access every element use getval to avoid creating a very inefficient full matrix!
objref m
m = new Matrix(3,4)
for i=0,m.nrow-1 {
for j=0, m.ncol-1 {
m.x[i][j] = 10*i + j
print i, j, m.x[i][j]
}
}
m.printf
xpanel("m")
xvalue("m(1,3) interpret", "m.x[1][3]", 1, "m.printf")
xpvalue("m(1,3) address", &m.x[1][3], 1, "m.printf")
xpanel()
Matrix
n = m.nrow
Matrixn = m.ncol
Matrix
mobj = msrcdest(nrow, ncol)objref m m = new Matrix(3,5) m for i=0,4 m.setcol(i,i) m.printf m.resize(7,7) m.printf() m.resize(4,2) m.printf()
Matrix
mdest = msrc.c()
Matrix
mdest = msrc.bcopy(i0, j0, n, m [, mout])mdest = msrc.bcopy(i0, j0, n, m, i1, j1 [, mout])objref m m = new Matrix(4,6) for i=0,m.nrow-1 for j=0,m.ncol-1 m.x[i][j] = 1 + 10*i + j m.printf m.bcopy(1,2,2,3).printf m.bcopy(1,2,2,3,2,3).printf m.bcopy(1,2,2,3,2,3, new Matrix(8,8)).printf
Matrix
val = m.getval(irow, jcol)
Matrix
0 = m.printf0 = m.printf("element_format")0 = m.printf("element_format", "row_format")
Matrix
0 = m.fprint(fileobj)0 = m.fprint(fileobj, "element_format")0 = m.fprint(fileobj, "element_format", "row_format")
Matrix
0 = m.scanf(File_object)0 = m.scanf(File_object, nrow, ncol)
objref m, f
f = new File("filename")
f.ropen()
m = new Matrix()
m.scanf(f)
print m.nrow, m.ncol
Matrix
vobj = msrc.mulv(vin)vobj = msrc.mulv(vin, vout)A sparse example execute following exampleprint "v1", v1 v1.printf print "m", m m.printf print "m*v1" m.mulv(v1).printf
objref m, v1
v1 = new Vector(100)
v1.indgen(1,1)
m = new Matrix(100, 100, 2) // sparse matrix
// reverse permutation
for i=0, 99 {
m.x[i][99 - i] = 1
}
m.mulv(v1).printf
Matrix
vobj = msrc.getrow(i)vobj = msrc.getrow(i, vout)
Matrix
vobj = msrc.getcol(i)vobj = msrc.getcol(i, vout)
Matrix
vobj = msrc.getdiag(i)vobj = msrc.getdiag(i, vout)
objref m
m = new Matrix(4,5)
for i=0, m.nrow-1 for j=0, m.ncol-1 m.x[i][j] = 1 + 10*j + 100*i
m.printf
for i=-m.nrow+1, m.ncol-1 {
printf("diagonal %d: ", i)
m.getdiag(i).printf
}
Matrix
vx = msrc.solv(vb)vx = msrc.solv(vb, vout and/or 1 in either order)Note: if the LUfactor is used, changes to the actual values of msrc would not affect the solution on subsequent calls to solv.
execute following exampleobjref m, b b = new Vector(3) b.indgen(1,1) m = new Matrix(3, 3) for i=0, m.nrow-1 for j=0, m.ncol-1 m.x[i][j] = i*j + 1 print "b" b.printf print "m" m.printf print "solution of m*x = b" m.solv(b).printf
objref m, b, x
m = new Matrix(1000, 1000, 2) // sparse type
m.setdiag(0, 3)
m.setdiag(-1, -1)
m.setdiag(1, -1)
b = new Vector(1000)
b.x[500] = 1
x = m.solv(b)
x.printf("%8.3f", 475, 525)
b.x[500] = 0
b.x[499] = 1
m.solv(b,1).printf("%8.3f", 475, 535)
Matrix
mobj = msrc.mulm(m)mobj = msrc.mulm(m, mout)
objref m1, m2, v1
m1 = new Matrix(6, 6)
for i=-1,1 {
if (i == 0) {
m1.setdiag(i, 2)
}else{
m1.setdiag(i, -1)
}
}
m2 = m1.inverse()
print "m1"
m1.printf
print "m2"
m2.printf(" %8.5f")
print "m1*m2"
m1.mulm(m2).printf(" %8.5f")
Matrix
mobj = m1srcdest.add(m2src)
Matrix
mobj = msrcdest.muls(scalar)objref m m = new Matrix(4,4) m.ident() m.muls(-10) m.printf
Matrix
mobj = msrcdest.setrow(i, vin)mobj = msrcdest.setrow(i, scalar)Otherwise fill the matrix row with a constant.
Matrix
mobj = msrcdest.setcol(i, vin)mobj = msrcdest.setcol(i, scalar)Otherwise fill the matrix column with a constant.
Matrix
mobj = msrcdest.setdiag(i, vin)mobj = msrcdest.setdiag(i, scalar)Otherwise fill the matrix diagonal with a constant.
objref v1, m
m = new Matrix(5,7)
v1 = new Vector(5)
for i=-4,6 {
m.setdiag(i, i)
}
m.printf
for i=-4,6 {
v1.indgen(1,1)
m.setdiag(i, v1)
}
m.printf
Matrix
mobj = msrcdest.zero()
Matrix
mobj = msrcdest.ident()objref m m = new Matrix(4,6) m.ident() m.printf()
Matrix
mobj = msrc.exp()mobj = msrc.exp(mout)
objref m, v1
m = new Matrix(8,8)
v1 = new Vector(8)
for i=-1,1 { v1.fill(2 - 3*abs(i)) m.setdiag(i, v1) }
m.exp().printf
Matrix
mobj = msrc.pow(i)mobj = msrc.pow(i, mout)
objref m
m = new Matrix(6, 6)
m.ident
m.x[0][5] = m.x[5][0] = 1
for i=0, 5 {
print i
m.pow(i).printf
}
Matrix
mobj = msrc.inverse()mobj = msrc.inverse(mout)
objref m, v1, minv
m = new Matrix(7,7)
v1 = new Vector(7)
for i=-1,1 { v1.fill(2 - 3*abs(i)) m.setdiag(i, v1) }
minv = m.inverse()
m.printf
minv.printf
m.mulm(minv).printf
Matrix
dvec = msrc.svd()dvec = msrc.svd(umat, vmat)
objref a, umat, vmat, dvec, dmat
proc svdtest() {
umat = new Matrix()
vmat = new Matrix()
dvec = $o1.svd(umat, vmat)
dmat = new Matrix($o1.nrow, $o1.ncol)
dmat.setdiag(0, dvec)
print "dvec" dvec.printf
print "dmat" dmat.printf
print "umat" umat.printf
print "vmat" vmat.printf
print "input ", $o1 $o1.printf()
print "ut*d*v"
umat.transpose.mulm(dmat).mulm(vmat).printf
}
a = new Matrix(5, 3)
a.setdiag(0, a.getdiag(0).indgen.add(1))
svdtest(a)
a = new Matrix(6, 6)
objref r
r = new Random()
r.discunif(1,10)
for i=0, a.nrow-1 {
a.setrow(i, a.getrow(i).setrand(r))
}
svdtest(a)
a = new Matrix(2,2)
a.setrow(0, 1)
a.setrow(1, 2)
svdtest(a)
Matrix
mdest = msrc.transpose()objref m m = new Matrix(1,5) for i=0, 4 m.x[0][i] = i m.printf m.transpose.printf m.transpose.mulm(m).printf m.mulm(m.transpose).printf
Matrix
veigenvalues = msrc.symmeig(eigenvectors)objref m, q, e m = new Matrix(5,5) m.setdiag(0, 2) m.setdiag(-1, -1) m.setdiag(1, -1) m.printf q = new Matrix(1,1) e = m.symmeig(q) print "eigenvectors" q.printf print "eigenvalues" e.printf print "qt*m*q" q.transpose.mulm(m).mulm(q).printf print "qt*q" q.transpose.mulm(q).printf
msrc must be symmetric but that fact is not checked.