| Class | AppMath::Mat |
| In: |
linalg.rb
|
| Parent: | Object |
Matrix space of arbitrary dimension. The intended usage is that the elements of a matrix are all either real or complex. Since one is allowed to change any matrix element into any object there is no guaranty for type-uniformity of the elements of a matrix.
| x | [RW] |
These are the 5 mehods to generate a matrix via ‘new‘
m1 = Mat.new(aMat) m2 = Mat.new(anArrayOfVec) m3 = Mat.new(aVec) m4 = Mat.new(aPositiveInteger, aRealOrComplex) m5 = Mat.new(aPositiveInteger, aPositiveInteger, aRealOrComplex)
Here, m1 is a copy of aMat, m2 is a matrix which has as row vectors, the components of anArrayOfVec. If these vectors have not all he same dimension, failure results; m3 is a square matrix in which only the main diagonal may have non-zero elements, and in which ths diagonal is given as aVec; m4 is a square matrix with the dimension given by the first argument, and with all matrix elements equal to the second argment; m5 is a rectangular matrix with dim1 and dim2 given by the first and the second argument, and with all matrix elements equal to the third argument.
# File linalg.rb, line 411
411: def initialize(*arg)
412: case arg.size
413: when 0
414: @x = Array.new
415: when 1 # ok if this is a Matrix or an array of Vectors
416: a0 = arg[0]
417: if a0.is_a?(Mat)
418: @x = Array.new(a0.x)
419: elsif a0.is_a?(Array)
420: n = a0.size
421: if n.zero?
422: @x = Array.new
423: else
424: misfit = 0
425: a0.each{|c|
426: misfit += 1 unless c.is_a?(Vec)
427: }
428: fail "input must consist of Vec-objects" unless misfit.zero?
429: misfit2 = 0
430: d2 = a0[1].dim
431: a0.each{|c|
432: misfit2 += 1 unless c.dim == d2
433: }
434: fail "input Vec-objects must agree in dimension" unless
435: misfit.zero?
436: @x = Array.new
437: a0.each{|c|
438: @x << c
439: }
440: end
441: elsif a0.is_a?(Vec) # make a diagonal matrix
442: n = a0.dim
443: if n.zero?
444: @x = Array.new
445: else
446: c = a0[1].to_0
447: vc = Vec.new(n,c)
448: @x = Array.new(n,vc)
449: for i in 1..n
450: s!(i,i,a0[i])
451: end
452:
453: end
454: else
455: fail "no reasonable construction available for this argument"
456: end
457: when 2 # make a square matrix, the diagonal filled with one element
458: # (all others zero)
459: n = arg[0]
460: a = arg[1]
461: zero = a.to_0
462: vc = Vec.new(n,zero)
463: @x = Array.new(n,vc)
464: for i in 1..n
465: vi = Vec.new(vc)
466: vi[i] = a
467: @x[i-1] = vi
468: end
469: when 3 # make rectangular matrix filled with one element
470: n1 = arg[0]
471: fail "first argument must be integer" unless n1.integer?
472: n2 = arg[1]
473: fail "second argument must be integer" unless n2.integer?
474: a = arg[2]
475: vc = Vec.new(n2,a)
476: @x = Array.new(n1,vc)
477: else
478: fail "no construction for more than 3 arguments"
479: end
480: end
Singular value decomposition. Slightly modified fom Press et al. Only needed as a algorithmic tool. The method for the end-user is method pseudo_inv.
# File linalg.rb, line 524
524: def Mat.svdcmp(a, w, v)
525:
526: m = a.dim1; n = a.dim2
527: fail "svdcmp: bad frame of a" unless m >= n
528: fail "svdcmp: bad frame of w" unless n == w.dim
529: fail "svdcmp: bad frame of v" unless v.dim1 == n && v.dim2 == n
530: fail "svdcmp: dim = 0 as input" if m.zero? || n.zero?
531:
532: iter_max=40
533: a11 = a[1][1]
534: zero =a11.to_0
535: one = a11.to_1
536: two = one + one
537: rv1 = Vec.new(n,zero)
538: g = zero; scale = zero; anorm = zero
539: for i in 1..n
540: l = i + 1
541: rv1[i] = scale * g
542: g = zero; s = zero; scale = zero
543: if i <= m
544: for k in i..m; scale += a[k][i].abs; end
545: if scale.nonzero?
546: for k in i..m
547: aki = a[k][i]
548: aki /= scale
549: a.s!(k,i,aki)
550: s += aki * aki
551: end
552: f = a[i][i]
553: g = - Basics.sign2(s.sqrt,f)
554: h = f * g - s
555: a.s!(i,i,f - g)
556: for j in l..n
557: s = zero
558: for k in i..m; s += a[k][i] * a[k][j]; end
559: f = s/h
560: for k in i..m
561: akj = a[k][j]
562: akj += f * a[k][i]
563: a.s!(k,j,akj)
564: end
565: end # for j in l..n
566: for k in i..m
567: aki = a[k][i]
568: aki *= scale
569: a.s!(k,i,aki)
570: end
571: end # scale != zero
572: end # i <= m
573: w[i] = scale * g
574: g = zero; s = R.c0; scale = zero
575: if i <= m && i != n
576: for k in l..n; scale += a[i][k].abs; end
577: if scale.nonzero?
578: for k in l..n
579: aik = a[i][k]
580: aik /= scale
581: a.s!(i,k,aik)
582: s += aik * aik
583: end
584: f = a[i][l]
585: g = - Basics.sign2(s.sqrt,f)
586: h = f * g - s
587: a.s!(i,l,f - g)
588: for k in l..n; rv1[k] = a[i][k]/h ; end
589: for j in l..m
590: s = zero
591: for k in l..n; s += a[j][k] * a[i][k]; end
592: for k in l..n
593: ajk = a[j][k]
594: ajk += s * rv1[k]
595: a.s!(j,k,ajk)
596: end
597: end # for j in l..m
598: for k in l..n; aik = a[i][k]; aik *= scale; a.s!(i,k,aik); end
599: end # if scale != zero
600: end # if i <= m && i != n
601: anorm = Basics.sup(anorm,w[i].abs + rv1[i].abs)
602: end # for i in 1..n
603: i = n
604: while i >= 1
605: if i < n
606: if g.nonzero?
607: for j in l..n; v.s!(j,i, (a[i][j]/a[i][l])/g); end
608: for j in l..n
609: s = zero
610: for k in l..n; s += a[i][k] * v[k][j]; end
611: for k in l..n
612: vkj =v[k][j]
613: vkj += s * v[k][i]
614: v.s!(k,j,vkj)
615: end
616: end # for j in l..n
617: end # if g.notzero!
618: for j in l..n; v.s!(i,j,zero); v.s!(j,i,zero); end
619: end # if i< n
620: v.s!(i,i,one)
621: g = rv1[i]
622: l = i
623: i -= 1
624: end # while i >= 1
625:
626: i = Basics.inf(m,n)
627: while i >= 1
628: l = i + 1
629: g = w[i]
630: for j in l..n; a.s!(i,j,zero); end
631: if g.nonzero?
632: g = one/g
633: for j in l..n
634: s = zero
635: for k in l..m; s += a[k][i] * a[k][j]; end
636: f = (s/a[i][i]) * g
637: for k in i..m
638: akj = a[k][j]; akj += f * a[k][i]; a.s!(k,j,akj)
639: end
640: end # for j in l..n
641: for j in i..m; aji = a[j][i]; aji *= g; a.s!(j,i,aji); end
642: else
643: for j in i..m; a.s!(j,i,zero); end
644: end # if g.nonzero?
645: aii = a[i][i]; aii += one; a.s!(i,i,aii)
646: i -= 1
647: end # while i >= 1
648:
649: k = n
650: while k >=1
651: for its in 1..iter_max
652: flag = 1
653: l = k
654: while l >= 1
655: nm = l - 1
656: if rv1[l].abs + anorm == anorm
657: flag = 0
658: break
659: end # if rv1[l].abs + anorm == anorm
660: break if w[nm].abs + anorm == anorm
661: l -= 1
662: end # while l >= 1
663: if flag.nonzero?
664: c = zero
665: s = one
666: for i in l..k
667: f = s * rv1[i]
668: rv1[i] = c * rv1[i]
669: break if f.abs + anorm == anorm
670: g = w[i]
671: h = f.hypot(g)
672: w[i] = h
673: h = one/h
674: c = g * h
675: s = -f * h
676: for j in 1..m
677: y = a[j][nm]; z = a[j][i];
678: a.s!(j,nm,y*c+z*s); a.s!(j,i,z*c-y*s)
679: end # for j in 1..m
680: end # for i in l..k
681: end # if flag.nonzero?
682: z = w[k]
683: if l == k
684: if z < zero
685: w[k] = -z
686: for j in 1..n; v.s!(j,k,-v[j][k]); end
687: end # if z < zero
688: break
689: end # if l == k
690: fail "no convergence in #{iter_max} iterations" if its == iter_max
691: x = w[l]; nm = k - 1; y = w[nm]; g = rv1[nm]; h = rv1[k]
692: f = ((y-z) * (y+z) + (g-h) * (g+h))/(h*y*two)
693: g = f.hypot(one)
694: f = ((x-z)*(x+z)+h*((y/(f+Basics.sign2(g,f)))-h))/x;
695: c = one; s = one
696: for j in l..nm
697: i = j + 1; g = rv1[i]; y =w[i]; h = s * g; g = c * g
698: z = f.hypot(h)
699: rv1[j] = z
700: c = f/z
701: s = h/z
702: f = x*c+g*s; g = g*c-x*s; h=y*s; y *= c;
703: for jj in 1..n
704: x=v[jj][j]; z=v[jj][i];
705: v.s!(jj,j,x*c+z*s); v.s!(jj,i,z*c-x*s)
706: end # for jj in 1..n
707: z = f.hypot(h)
708: w[j] = z
709: if z.nonzero?
710: z = one/z; c = f * z; s = h * z
711: end # if z.nonzero?
712: f=c*g+s*y; x=c*y-s*g;
713: for jj in 1..m
714: y=a[jj][j]; z=a[jj][i]
715: a.s!(jj,j,y*c+z*s); a.s!(jj,i,z*c-y*s)
716: end # for jj in 1..m
717: end # for j in l..nm
718: rv1[l] = zero
719: rv1[k] = f
720: w[k] = x
721: end # for its in 1..iter_max
722: k -= 1
723: end # while k >=1
724: end
Consistency test of class Mat.
# File linalg.rb, line 1038
1038: def Mat.test(n0, verbose = true, complex = false )
1039: puts "Doing Mat.test( n = #{n0}, verbose = #{verbose}," +
1040: " complex = #{complex} ) for R.prec = #{R.prec}:"
1041: puts "******************************************************"
1042: t1 = Time.now
1043: s = R.c0
1044: puts "class of s is " + s.class.to_s
1045: i = n0 + 137
1046: a = Mat.tob(n0, i, complex)
1047: i += 1
1048: b = Mat.tob(n0, i, complex)
1049: i += 1
1050: c = Mat.tob(n0, i, complex)
1051: x = Vec.tob(n0, i, complex)
1052: i += 1
1053: y = Vec.tob(n0, i, complex)
1054: i += 1
1055: s1 = complex ? C.ran(i) : R.ran(i)
1056: i += 1
1057: s2 = complex ? C.ran(i) : R.ran(i)
1058:
1059: unit = a.to_1
1060:
1061: a0 = a.clone
1062: b0 = b.clone
1063: c0 = c.clone
1064:
1065: abc0 = a0.abs + b0.abs + c0.abs
1066:
1067: ac = a.clone
1068: a = b
1069: r = a
1070: l = b
1071: ds = r.dis(l)
1072: puts "assignment of variables: ds = " + ds.to_s if verbose
1073: s += ds
1074:
1075: a = ac
1076: r = a
1077: l = ac
1078: ds = r.dis(l)
1079: puts "assignment of variables 2: ds = " + ds.to_s if verbose
1080: s += ds
1081:
1082: r = (a + b) + c
1083: l = a + (b + c)
1084: ds = r.dis(l)
1085: puts "associativity of +: ds = " + ds.to_s if verbose
1086: s += ds
1087:
1088: r = (a - b) + c
1089: l = a - (b - c)
1090: ds = r.dis(l)
1091: puts "associativity of -: ds = " + ds.to_s if verbose
1092: s += ds
1093:
1094: r = (a * b) * c
1095: l = a * (b * c)
1096: ds = r.dis(l)
1097: puts "associativity of *: ds = " + ds.to_s if verbose
1098: s += ds
1099:
1100: r = (a + b) * s1
1101: l = a * s1 + b * s1
1102: ds = r.dis(l)
1103: puts "distributivity of multiplication by scalars: ds = " +
1104: ds.to_s if verbose
1105: s += ds
1106:
1107: r = c * (s1*s2)
1108: l = (c * s1) * s2
1109: ds = r.dis(l)
1110: puts "distributivity of multiplication by scalars: ds = " +
1111: ds.to_s if verbose
1112: s += ds
1113:
1114: r = a
1115: l = -(-a)
1116: ds = r.dis(l)
1117: puts "idempotency of unary minus: ds = " + ds.to_s if verbose
1118: s += ds
1119:
1120: r = (a + b).spr(c)
1121: l = a.spr(c) + b.spr(c)
1122: ds = r.dis(l)
1123: puts "distributivity of spr: ds = " + ds.to_s if verbose
1124: s += ds
1125:
1126: r = (a + b) * c
1127: l = a * c + b * c
1128: ds = r.dis(l)
1129: puts "distributivity of matrix multiplication: ds = " +
1130: ds.to_s if verbose
1131: s += ds
1132:
1133: r = (a * b) * x
1134: l = a * (b * x)
1135: ds = r.dis(l)
1136: puts "action on vectors 1: ds = " + ds.to_s if verbose
1137: s += ds
1138:
1139: r = (a + b) * x
1140: l = a * x + b * x
1141: ds = r.dis(l)
1142: puts "action on vectors 2: ds = " + ds.to_s if verbose
1143: s += ds
1144:
1145: r = b * (x + y)
1146: l = b * x + b * y
1147: ds = r.dis(l)
1148: puts "action on vectors 3: ds = " + ds.to_s if verbose
1149: s += ds
1150:
1151: r = c * (x * s1)
1152: l = (c * s1) * x
1153: ds = r.dis(l)
1154: puts "action on vectors 4: ds = " + ds.to_s if verbose
1155: s += ds
1156:
1157: if complex == false
1158: r = unit
1159: l = a * a.pseudo_inv
1160: ds = r.dis(l)
1161: puts "pseudo inverse is right inverse: ds = " + ds.to_s if verbose
1162: s += ds
1163:
1164: r = unit
1165: l = a.pseudo_inv * a
1166: ds = r.dis(l)
1167: puts "pseudo inverse is left inverse: ds = " + ds.to_s if verbose
1168: s += ds
1169: else
1170: puts "test of pseudo inverse left out, since not implemented for complex"
1171: end
1172:
1173: aMem = a.clone
1174: bMem = b.clone
1175: cMem = c.clone
1176:
1177: # Testing the access functions, under harsh conditions with
1178: # inserting double transposition
1179:
1180: for i in 1..a.dim1
1181: for j in 1..a.dim2
1182: b.s!(i,j,a[i][j]) # copies a to b
1183: end
1184: end
1185:
1186: l = a
1187: r = a.trp.trp
1188: ds = r.dis(l)
1189: puts "test of double transposition: ds = " + ds.to_s if verbose
1190: s += ds
1191:
1192: for i in 1..b.dim1
1193: for j in 1..b.dim2
1194: c.s!(i,j,b[i][j]) # copies b to c
1195: end
1196: end
1197: # Finally c should have the content of a
1198:
1199: r = a
1200: l = c
1201: ds = r.dis(l)
1202: puts "test of access functions: ds = " + ds.to_s if verbose
1203: s += ds
1204:
1205: a = aMem; b = bMem; c = cMem
1206:
1207: abc = a.abs + b.abs + c.abs
1208:
1209: l = abc
1210: r = abc0
1211: ds = r.dis(l)
1212: puts "heavy test of assignment: ds = " + ds.to_s if verbose
1213:
1214: t2 = Time.now
1215:
1216: if verbose
1217: puts
1218: a.prn("a")
1219: puts
1220: b.prn("b")
1221: puts
1222: c.prn("c")
1223: puts
1224: s1.prn("s1")
1225: puts
1226: s2.prn("s2")
1227: puts
1228: x.prn("x")
1229: puts
1230: y.prn("y")
1231: end
1232:
1233: puts "class of s is " + s.class.to_s + " ."
1234: puts "The error sum s is " + s.to_s + " ."
1235: puts "It should be close to 0."
1236: puts "Computation time was " + (t2-t1).to_s
1237: s
1238: end
Generates a test object, here a n times n matrix with random elements. This object depends rather chaotically on the integer parameter i. If the last argument is ‘false’ the test matrix will have R-typed elements, and C-typed elements else.
# File linalg.rb, line 492
492: def Mat.tob(n,i, complex = false)
493: if complex
494: ri = C.tob(i)
495: zero = ri.to_0
496: res=Mat.new(n, n, zero)
497: rg1 = Ran.new(-ri.re,ri.re)
498: rg2 = Ran.new(-ri.im,ri.im)
499: for j in 1..n
500: for k in 1..n
501: yjk = C.new(rg1.ran,rg2.ran)
502: res.s!(j,k,yjk)
503: end
504: end
505: res
506: else
507: ri = R.tob(i)
508: zero = ri.to_0
509: res=Mat.new(n, n, zero)
510: rg = Ran.new(-ri,ri)
511: for j in 1..n
512: for k in 1..n
513: yjk = rg.ran
514: res.s!(j,k,yjk)
515: end
516: end
517: res
518: end
519: end
Multiplication of a Mat with either a Mat, Vec, or Numeric
# File linalg.rb, line 861
861: def *(v)
862: d1 = dim1
863: fail "dim1 == 0" if d1.zero?
864: d2 = dim2
865: fail "dim2 == 0" if d2.zero?
866: zero=@x[0].x[0].to_0
867: if v.is_a?(Mat)
868: d3 = v.dim1
869: fail "dimenson mismatch" unless d3 == d2
870: d4 = v.dim2
871: fail "dim4 == 0" if d4.zero?
872: res = Mat.new(d1,d4,zero)
873: for i in 0...d1
874: ip = i + 1
875: for j in 0...d4
876: vij = zero
877: jp = j + 1
878: for k in 0...d2
879: vij += @x[i].x[k] * v.x[k].x[j]
880: end
881: res.s!(ip, jp, vij)
882: end
883: end
884: elsif v.is_a?(Vec)
885: d3 = v.dim
886: fail "dimenson mismatch" unless d3 == d2
887: res = Vec.new(d1,zero)
888: for i in 0...d1
889: vi = zero
890: for j in 0...d2
891: vi += @x[i].x[j] * v.x[j]
892: end
893: res.x[i] = vi
894: end
895: elsif v.is_a?(Numeric) # multiplication with scalar
896: res = clone
897: for i in 1..d1
898: res[i] *= v
899: end
900: else
901: fail "can't multiply with this object"
902: end
903: res
904: end
Unary minus operator. Returns - self.
# File linalg.rb, line 770
770: def -@
771: res = clone
772: for i in 1..dim
773: res[i] = -res[i]
774: end
775: res
776: end
The order relation is here lexicographic ordering of lists. Needed only for book-keeping purposes.
# File linalg.rb, line 908
908: def <=> (v)
909: d1 = dim; d2 = v.dim
910: if d1 < d2
911: return -1
912: elsif d1 > d2
913: return 1
914: else
915: for i in 1..d1
916: ci = self[i] <=> v[i]
917: return ci unless ci == 0
918: end
919: end
920: return 0
921: end
Absolute value, always real
# File linalg.rb, line 944
944: def abs
945: if complex?
946: abs2.re.sqrt
947: else
948: abs2.sqrt
949: end
950: end
Returns an independent copy of self.
# File linalg.rb, line 483
483: def clone
484: Mat.new(self)
485: end
# File linalg.rb, line 967
967: def complex?
968: return nil if dim.zero?
969: @x[0].complex?
970: end
Returns the Hermitian conjugate of self.
# File linalg.rb, line 818
818: def conj
819: d1 = dim1
820: d2 = dim2
821: fail "dim1 == 0" if d1.zero?
822: fail "dim2 == 0" if d2.zero?
823: zero = @x[0].x[0].to_0
824: res = Mat.new(d2,d1,zero)
825: for i in 0...d1
826: for j in 0...d2
827: sij = @x[i].x[j]
828: res.s!(j+1,i+1,sij.conj)
829: end
830: end
831: res
832: end
Returns the ‘dimension’ of the matrix, i.e. the number of its row-vectors. This thus is m for a ‘m times n - matrix’.
# File linalg.rb, line 382
382: def dim; @x.size; end
Relative distance between matrices
# File linalg.rb, line 953
953: def dis(v)
954: a = abs
955: b = v.abs
956: d = (self - v).abs
957: s = a + b
958: return R.c0 if s.zero?
959: d1 = d/s
960: d < d1 ? d : d1
961: end
In most cases (‘up to a subst of measure zero’) the pseudo-inverse is also the inverse.
# File linalg.rb, line 1033
1033: def inv
1034: pseudo_inv
1035: end
Prints the content of self and naming the output.
# File linalg.rb, line 761
761: def prn(name)
762: for i in 1..dim1
763: for j in 1..dim2
764: puts " #{name}[#{i}][#{j}] = " + self[i][j].to_s
765: end
766: end
767: end
Returns the pseudo-inverse (als known as Penrose inverse) of self. If the argument acc is not zero, the discontinous treatment of singular values near zero is replaced by a continuous one. Notice that the pseudo inverse always exists, and that the pseudo- inverse of a (m,n)-matrix is a (n,m)-matrix. If the argument acc is not zero, the pseudo-inverse is the first component of a 4-array res, which also contains the the intermediary quantities a, w, v resulting from the call
Mat.svdcmp(a,w,v). a == res[1], w == res[2], v == res[3].
Especially the list w of the original singular values is thus made accessible, so that one can judge whether their processing controlled by the parameter acc was reasonable. The pseudo_inverse is the most useful and stable mehod to solve linear equations: Let a be a (m,n)-matrix, and b a m-vector. The equation
a * x = b (i)
determines a n-vector x as
x = a.pseudo_inv * b
which is a solution of (i) if there is one. If there are many solutions it is the one of minimum absolute value, and if there is no solution it comes closest to be a solution: It minimizes the defect
(a * x - b).abs
Simply great! No need for LU-decompositions any more.
# File linalg.rb, line 996
996: def pseudo_inv(acc = 0)
997: if complex?
998: fail "Pseudo inverse not yet impemented for complex matrices"
999: end
1000: m = dim1
1001: fail "dim1 == 0" if m.zero?
1002: n=dim2
1003: fail "dim2 == 0" if n.zero?
1004: rightframe = m >= n
1005: a = clone
1006: a = a.trp if rightframe == false
1007: m = a.dim1; n = a.dim2
1008: zero = a[1][1].to_0
1009: v = Mat.new(n,n,zero)
1010: w = Vec.new(n,zero)
1011: vr = Mat.new(n,m,zero) # sic
1012: Mat.svdcmp(a,w,v)
1013: wi = w.pseudo_inv(acc) # w does not come out orderd
1014: for i in 1..n
1015: for j in 1..m
1016: sum = zero
1017: for k in 1..n
1018: sum += v[i][k] * wi[k] * a[j][k]
1019: end
1020: vr.s!(i,j,sum)
1021: end
1022: end
1023: vr = vr.trp if rightframe == false
1024: if acc.zero?
1025: vr
1026: else
1027: [vr,a,w,v]
1028: end
1029: end
The ‘s’ stands for ‘set (value)’ and the ’!’ this is a method by which self changes (non-constant or mutating method).
# File linalg.rb, line 741
741: def s!(i,j,a)
742: si = Vec.new(self[i]) # don't change the row-vector self[i]
743: # itself. Such changes are subject to subtle side effects.
744: # Worcing on a copy is safe.
745: si[j] = a # changing si here is normal syntax
746: self[i] = si # this is OK, of course
747: # self[i] = Vec.new(si) would work too, but would cause more work
748: end
Scalar product of matrices.
# File linalg.rb, line 928
928: def spr(v)
929: fail "dimension mismatch" unless dim == v.dim
930: return nil if dim.zero?
931: s = self[1].spr(v[1])
932: for i in 2..dim
933: s += self[i].spr(v[i])
934: end
935: s
936: end
‘to zero’ Returns a matrix with he same dimensions as self, but with all matrix elements set to zero.
# File linalg.rb, line 838
838: def to_0
839: d1 = dim1
840: d2 = dim2
841: fail "dim1 == 0" if d1.zero?
842: fail "dim2 == 0" if d2.zero?
843: zero = @x[0].x[0].to_0
844: Mat.new(d2,d1,zero)
845: end
‘to one’. Returns a matrix with he same dimensions as self, but with all matrix elements set to 1.
# File linalg.rb, line 849
849: def to_1
850: d1 = dim1
851: d2 = dim2
852: fail "dim1 == 0" if d1.zero?
853: fail "dim2 == 0" if d2.zero?
854: fail "dim1 != dim2" unless d1 == d2
855: unit = @x[0].x[0].to_1
856: diag = Vec.new(d1,unit)
857: Mat.new(diag)
858: end
Returns a string which consists of a list of the strings which represent the row-vectors.
# File linalg.rb, line 752
752: def to_s
753: res = "\n Mat"
754: for i in 0...dim
755: res += "\n " + x[i].to_s
756: end
757: res + "\n end Mat"
758: end
Returns the transposed of self.
# File linalg.rb, line 801
801: def trp
802: d1 = dim1
803: d2 = dim2
804: fail "dim1 == 0" if d1.zero?
805: fail "dim2 == 0" if d2.zero?
806: zero = @x[0].x[0].to_0
807: res = Mat.new(d2,d1,zero)
808: for i in 0...d1
809: for j in 0...d2
810: sij = @x[i].x[j]
811: res.s!(j+1,i+1,sij)
812: end
813: end
814: res
815: end