classes add getdiag nrow solv bcopy getrow pow spgetrowval c getval printf sprowlen cholesky_factor ident resize svd det inverse scanf symmeig exp mulm setcol to_vector fprint muls setdiag transpose from_vector mulv setrow x getcol ncol setval zero

`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] = val`

`expr(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()

In Python, the m.x[i][j] syntax does not work and one must use the setval function

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

`val = m.setval(irow, jcol, val)`

Matrix

`n = m.sprowlen(i)`

Matrix

`x = m.spgetrowval(i, jx, &j)`

proc sparse_print() { local i, j, jx, x print $o1 for i=0, $o1.nrow-1 { printf("%d ", i) for jx = 0, $o1.sprowlen(i)-1 { x = $o1.spgetrowval(i, jx, &j) printf(" %d:%g", j, x) } printf ("\n") } } objref m m = new Matrix(4, 5, 2) m.x[0][2] = 1.2 m.x[0][4] = 2.4 m.x[1][1] = 3.1 for i=0, 4 { m.x[3][i] = i/10 } sparse_print(m)

Matrix

`0 = m.printf`

`0 = 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")`

`0 = m.fprint(0, fileobj [,...])`

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

`mantissa = m.det(&base10exponent)`

objref m m = new Matrix(2,2) m.x[0][1] = 20 m.x[1][0] = 30 m.printf() ex = 0 mant = m.det(&ex) print mant*10^ex

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.

Matrix

`vobj = msrc.to_vector()`

`vobj = msrc.to_vector(vout)`

objref m m = new Matrix(4,5) m.from_vector(m.to_vector().indgen).printf

Matrix

`mobj = msrcdest.from_vector(vec)`

objref m m = new Matrix(4,5) m.from_vector(m.to_vector().indgen).printf

Matrix

`mc = msrcdest.cholesky_factor()`

objref m, cf m = new Matrix(3,3) for i=0,2 for j=0,2 m.x[i][j] = (i+j)*(i+j) m.printf cf = m.c.cholesky_factor() cf.mulm(cf.transpose()).printf