diff --git a/src/numeric.js b/src/numeric.js index 537b68f..1a115c6 100644 --- a/src/numeric.js +++ b/src/numeric.js @@ -22,7 +22,7 @@ numeric.bench = function bench (f,interval) { while(i>0) { f(); i--; } t2 = new Date(); return 1000*(3*n-1)/(t2-t1); -} +}; numeric._myIndexOf = (function _myIndexOf(w) { var n = this.length,k; @@ -95,7 +95,7 @@ numeric.prettyPrint = function prettyPrint(x) { } foo(x); return ret.join(''); -} +}; numeric.parseDate = function parseDate(d) { function foo(d) { @@ -106,7 +106,7 @@ numeric.parseDate = function parseDate(d) { return ret; } return foo(d); -} +}; numeric.parseFloat = function parseFloat_(d) { function foo(d) { @@ -117,7 +117,7 @@ numeric.parseFloat = function parseFloat_(d) { return ret; } return foo(d); -} +}; numeric.parseCSV = function parseCSV(t) { var foo = t.split('\n'); @@ -125,7 +125,7 @@ numeric.parseCSV = function parseCSV(t) { var ret = []; var pat = /(([^'",]*)|('[^']*')|("[^"]*")),/g; var patnum = /^\s*(([+-]?[0-9]+(\.[0-9]*)?(e[+-]?[0-9]+)?)|([+-]?[0-9]*(\.[0-9]+)?(e[+-]?[0-9]+)?))\s*$/; - var stripper = function(n) { return n.substr(0,n.length-1); } + var stripper = function(n) { return n.substr(0,n.length-1); }; var count = 0; for(k=0;k=0;i--) { ret[i] = numeric.rep(s,v,k+1); } return ret; -} +}; numeric.dotMMsmall = function dotMMsmall(x,y) { @@ -413,7 +413,7 @@ numeric.dotMMsmall = function dotMMsmall(x,y) { ret[i] = foo; } return ret; -} +}; numeric._getCol = function _getCol(A,j,x) { var n = A.length, i; for(i=n-1;i>0;--i) { @@ -422,7 +422,7 @@ numeric._getCol = function _getCol(A,j,x) { x[i] = A[i][j]; } if(i===0) x[0] = A[0][j]; -} +}; numeric.dotMMbig = function dotMMbig(x,y){ var gc = numeric._getCol, p = y.length, v = Array(p); var m = x.length, n = y[0].length, A = new Array(m), xj; @@ -441,14 +441,14 @@ numeric.dotMMbig = function dotMMbig(x,y){ } } return A; -} +}; numeric.dotMV = function dotMV(x,y) { var p = x.length, q = y.length,i; var ret = Array(p), dotVV = numeric.dotVV; for(i=p-1;i>=0;i--) { ret[i] = dotVV(x[i],y); } return ret; -} +}; numeric.dotVM = function dotVM(x,y) { var i,j,k,p,q,r,ret,foo,bar,woo,i0,k0,p0,r0,s1,s2,s3,baz,accum; @@ -464,7 +464,7 @@ numeric.dotVM = function dotVM(x,y) { ret[k] = woo; } return ret; -} +}; numeric.dotVV = function dotVV(x,y) { var i,n=x.length,i1,ret = x[n-1]*y[n-1]; @@ -474,7 +474,7 @@ numeric.dotVV = function dotVV(x,y) { } if(i===0) { ret += x[0]*y[0]; } return ret; -} +}; numeric.dot = function dot(x,y) { var d = numeric.dim; @@ -490,7 +490,7 @@ numeric.dot = function dot(x,y) { case 0: return x*y; default: throw new Error('numeric.dot only works on vectors and matrices'); } -} +}; numeric.diag = function diag(d) { var i,i1,j,n = d.length, A = Array(n), Ai; @@ -511,7 +511,7 @@ numeric.diag = function diag(d) { A[i] = Ai; } return A; -} +}; numeric.getDiag = function(A) { var n = Math.min(A.length,A[0].length),i,ret = Array(n); for(i=n-1;i>=1;--i) { @@ -523,9 +523,9 @@ numeric.getDiag = function(A) { ret[0] = A[0][0]; } return ret; -} +}; -numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); } +numeric.identity = function identity(n) { return numeric.diag(numeric.rep([n],1)); }; numeric.pointwise = function pointwise(params,body,setup) { if(typeof setup === "undefined") { setup = ""; } var fun = []; @@ -558,7 +558,7 @@ numeric.pointwise = function pointwise(params,body,setup) { 'return ret;' ); return Function.apply(null,fun); -} +}; numeric.pointwise2 = function pointwise2(params,body,setup) { if(typeof setup === "undefined") { setup = ""; } var fun = []; @@ -583,7 +583,7 @@ numeric.pointwise2 = function pointwise2(params,body,setup) { 'return ret;' ); return Function.apply(null,fun); -} +}; numeric._biforeach = (function _biforeach(x,y,s,k,f) { if(k === s.length-1) { f(x,y); return; } var i,n=s[k]; @@ -790,7 +790,7 @@ numeric.trunc = function trunc(x,y) { } if (typeof y === "object") return numeric.truncSV(x,y); return Math.round(x/y)*y; -} +}; numeric.inv = function inv(x) { var s = numeric.dim(x), abs = Math.abs, m = s[0], n = s[1]; @@ -818,7 +818,7 @@ numeric.inv = function inv(x) { } } return I; -} +}; numeric.det = function det(x) { var s = numeric.dim(x); @@ -846,7 +846,7 @@ numeric.det = function det(x) { ret *= Aj[j]; } return ret*A[j][j]; -} +}; numeric.transpose = function transpose(x) { var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj; @@ -873,7 +873,7 @@ numeric.transpose = function transpose(x) { if(j===0) { ret[0][0] = A0[0]; } } return ret; -} +}; numeric.negtranspose = function negtranspose(x) { var i,j,m = x.length,n = x[0].length, ret=Array(n),A0,A1,Bj; for(j=0;j=0;i--) ret[i] = _random(s,k+1); return ret; -} -numeric.random = function random(s) { return numeric._random(s,0); } +}; +numeric.random = function random(s) { return numeric._random(s,0); }; -numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); } +numeric.norm2 = function norm2(x) { return Math.sqrt(numeric.norm2Squared(x)); }; numeric.linspace = function linspace(a,b,n) { if(typeof n === "undefined") n = Math.max(Math.round(b-a)+1,1); @@ -926,7 +926,7 @@ numeric.linspace = function linspace(a,b,n) { n--; for(i=n;i>=0;i--) { ret[i] = (i*b+(n-i)*a)/n; } return ret; -} +}; numeric.getBlock = function getBlock(x,from,to) { var s = numeric.dim(x); @@ -940,7 +940,7 @@ numeric.getBlock = function getBlock(x,from,to) { return ret; } return foo(x,0); -} +}; numeric.setBlock = function setBlock(x,from,to,B) { var s = numeric.dim(x); @@ -951,7 +951,7 @@ numeric.setBlock = function setBlock(x,from,to,B) { } foo(x,B,0); return x; -} +}; numeric.getRange = function getRange(A,I,J) { var m = I.length, n = J.length; @@ -964,7 +964,7 @@ numeric.getRange = function getRange(A,I,J) { for(j=n-1;j!==-1;--j) Bi[j] = AI[J[j]]; } return B; -} +}; numeric.blockMatrix = function blockMatrix(X) { var s = numeric.dim(X); @@ -990,7 +990,7 @@ numeric.blockMatrix = function blockMatrix(X) { I += X[i][0].length; } return Z; -} +}; numeric.tensor = function tensor(x,y) { if(typeof x === "number" || typeof y === "number") return numeric.mul(x,y); @@ -1015,11 +1015,11 @@ numeric.tensor = function tensor(x,y) { A[i] = Ai; } return A; -} +}; // 3. The Tensor type T -numeric.T = function T(x,y) { this.x = x; this.y = y; } -numeric.t = function t(x,y) { return new numeric.T(x,y); } +numeric.T = function T(x,y) { this.x = x; this.y = y; }; +numeric.t = function t(x,y) { return new numeric.T(x,y); }; numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) { var io = numeric.indexOf; @@ -1047,7 +1047,7 @@ numeric.Tbinop = function Tbinop(rr,rc,cr,cc,setup) { '}\n'+ 'return new numeric.T('+rr+');\n' ); -} +}; numeric.T.prototype.add = numeric.Tbinop( 'add(x.x,y.x)', @@ -1072,14 +1072,14 @@ numeric.T.prototype.reciprocal = function reciprocal() { return new numeric.T(div(this.x,d),div(numeric.neg(this.y),d)); } return new T(div(1,this.x)); -} +}; numeric.T.prototype.div = function div(y) { if(!(y instanceof numeric.T)) y = new numeric.T(y); if(y.y) { return this.mul(y.reciprocal()); } var div = numeric.div; if(this.y) { return new numeric.T(div(this.x,y.x),div(this.y,y.x)); } return new numeric.T(div(this.x,y.x)); -} +}; numeric.T.prototype.dot = numeric.Tbinop( 'dot(x.x,y.x)', 'dot(x.x,y.x),dot(x.x,y.y)', @@ -1090,12 +1090,12 @@ numeric.T.prototype.transpose = function transpose() { var t = numeric.transpose, x = this.x, y = this.y; if(y) { return new numeric.T(t(x),t(y)); } return new numeric.T(t(x)); -} +}; numeric.T.prototype.transjugate = function transjugate() { var t = numeric.transpose, x = this.x, y = this.y; if(y) { return new numeric.T(t(x),numeric.negtranspose(y)); } return new numeric.T(t(x)); -} +}; numeric.Tunop = function Tunop(r,c,s) { if(typeof s !== "string") { s = ''; } return Function( @@ -1106,7 +1106,7 @@ numeric.Tunop = function Tunop(r,c,s) { '}\n'+ r+';\n' ); -} +}; numeric.T.prototype.exp = numeric.Tunop( 'return new numeric.T(ex)', @@ -1202,7 +1202,7 @@ numeric.T.prototype.inv = function inv() { } } return new numeric.T(Rx,Ry); -} +}; numeric.T.prototype.get = function get(i) { var x = this.x, y = this.y, k = 0, ik, n = i.length; if(y) { @@ -1220,7 +1220,7 @@ numeric.T.prototype.get = function get(i) { k++; } return new numeric.T(x); -} +}; numeric.T.prototype.set = function set(i,v) { var x = this.x, y = this.y, k = 0, ik, n = i.length, vx = v.x, vy = v.y; if(n===0) { @@ -1267,7 +1267,7 @@ numeric.T.prototype.set = function set(i,v) { ik = i[k]; x[ik] = vx; return this; -} +}; numeric.T.prototype.getRows = function getRows(i0,i1) { var n = i1-i0+1, j; var rx = Array(n), ry, x = this.x, y = this.y; @@ -1278,7 +1278,7 @@ numeric.T.prototype.getRows = function getRows(i0,i1) { return new numeric.T(rx,ry); } return new numeric.T(rx); -} +}; numeric.T.prototype.setRows = function setRows(i0,i1,A) { var j; var rx = this.x, ry = this.y, x = A.x, y = A.y; @@ -1290,12 +1290,12 @@ numeric.T.prototype.setRows = function setRows(i0,i1,A) { for(j=i0;j<=i1;j++) { ry[j] = numeric.rep([x[j-i0].length],0); } } return this; -} +}; numeric.T.prototype.getRow = function getRow(k) { var x = this.x, y = this.y; if(y) { return new numeric.T(x[k],y[k]); } return new numeric.T(x[k]); -} +}; numeric.T.prototype.setRow = function setRow(i,v) { var rx = this.x, ry = this.y, x = v.x, y = v.y; rx[i] = x; @@ -1306,13 +1306,13 @@ numeric.T.prototype.setRow = function setRow(i,v) { ry = numeric.rep([x.length],0); } return this; -} +}; numeric.T.prototype.getBlock = function getBlock(from,to) { var x = this.x, y = this.y, b = numeric.getBlock; if(y) { return new numeric.T(b(x,from,to),b(y,from,to)); } return new numeric.T(b(x,from,to)); -} +}; numeric.T.prototype.setBlock = function setBlock(from,to,A) { if(!(A instanceof numeric.T)) A = new numeric.T(A); var x = this.x, y = this.y, b = numeric.setBlock, Ax = A.x, Ay = A.y; @@ -1324,31 +1324,31 @@ numeric.T.prototype.setBlock = function setBlock(from,to,A) { } b(x,from,to,Ax); if(y) b(y,from,to,numeric.rep(numeric.dim(Ax),0)); -} +}; numeric.T.rep = function rep(s,v) { var T = numeric.T; if(!(v instanceof T)) v = new T(v); var x = v.x, y = v.y, r = numeric.rep; if(y) return new T(r(s,x),r(s,y)); return new T(r(s,x)); -} +}; numeric.T.diag = function diag(d) { if(!(d instanceof numeric.T)) d = new numeric.T(d); var x = d.x, y = d.y, diag = numeric.diag; if(y) return new numeric.T(diag(x),diag(y)); return new numeric.T(diag(x)); -} +}; numeric.T.eig = function eig() { if(this.y) { throw new Error('eig: not implemented for complex matrices.'); } return numeric.eig(this.x); -} -numeric.T.identity = function identity(n) { return new numeric.T(numeric.identity(n)); } +}; +numeric.T.identity = function identity(n) { return new numeric.T(numeric.identity(n)); }; numeric.T.prototype.getDiag = function getDiag() { var n = numeric; var x = this.x, y = this.y; if(y) { return new n.T(n.getDiag(x),n.getDiag(y)); } return new n.T(n.getDiag(x)); -} +}; // 4. Eigenvalues of real matrices @@ -1360,7 +1360,7 @@ numeric.house = function house(x) { var foo = numeric.norm2(v); if(foo === 0) { /* this should not happen */ throw new Error('eig: internal error'); } return numeric.div(v,foo); -} +}; numeric.toUpperHessenberg = function toUpperHessenberg(me) { var s = numeric.dim(me); @@ -1384,7 +1384,7 @@ numeric.toUpperHessenberg = function toUpperHessenberg(me) { } } return {H:A, Q:Q}; -} +}; numeric.epsilon = 2.220446049250313e-16; @@ -1473,7 +1473,7 @@ numeric.QRFrancis = function(H,maxiter) { } } throw new Error('numeric: eigenvalue iteration does not converge -- increase maxiter?'); -} +}; numeric.eig = function eig(A,maxiter) { var QH = numeric.toUpperHessenberg(A); @@ -1589,7 +1589,7 @@ numeric.ccsSparse = function ccsSparse(A) { } } return [Ai,Aj,Av]; -} +}; numeric.ccsFull = function ccsFull(A) { var Ai = A[0], Aj = A[1], Av = A[2], s = numeric.ccsDim(A), m = s[0], n = s[1], i,j,j0,j1,k; var B = numeric.rep([m,n],0); @@ -1599,7 +1599,7 @@ numeric.ccsFull = function ccsFull(A) { for(j=j0;jn) k.pop(); return ret; -} +}; // 6. Coordinate matrices numeric.cLU = function LU(A) { @@ -2338,7 +2338,7 @@ numeric.cgrid = function grid(n,shape) { if(typeof shape !== "function") { switch(shape) { case 'L': - shape = function(i,j) { return (i>=n[0]/2 || j=n[0]/2 || j b) - return a*Math.sqrt(1.0+(b*b/a/a)) - else if (b == 0.0) - return a - return b*Math.sqrt(1.0+(a*a/b/b)) - } - - //Householder's reduction to bidiagonal form - - var f= 0.0; - var g= 0.0; - var h= 0.0; - var x= 0.0; - var y= 0.0; - var z= 0.0; - var s= 0.0; - - for (i=0; i < n; i++) - { - e[i]= g; - s= 0.0; - l= i+1; - for (j=i; j < m; j++) - s += (u[j][i]*u[j][i]); - if (s <= tolerance) - g= 0.0; - else - { - f= u[i][i]; - g= Math.sqrt(s); - if (f >= 0.0) g= -g; - h= f*g-s - u[i][i]=f-g; - for (j=l; j < n; j++) - { - s= 0.0 - for (k=i; k < m; k++) - s += u[k][i]*u[k][j] - f= s/h - for (k=i; k < m; k++) - u[k][j]+=f*u[k][i] - } - } - q[i]= g - s= 0.0 - for (j=l; j < n; j++) - s= s + u[i][j]*u[i][j] - if (s <= tolerance) - g= 0.0 - else - { - f= u[i][i+1] - g= Math.sqrt(s) - if (f >= 0.0) g= -g - h= f*g - s - u[i][i+1] = f-g; - for (j=l; j < n; j++) e[j]= u[i][j]/h - for (j=l; j < m; j++) - { - s=0.0 - for (k=l; k < n; k++) - s += (u[j][k]*u[i][k]) - for (k=l; k < n; k++) - u[j][k]+=s*e[k] - } - } - y= Math.abs(q[i])+Math.abs(e[i]) - if (y>x) - x=y - } - - // accumulation of right hand gtransformations - for (i=n-1; i != -1; i+= -1) - { - if (g != 0.0) - { - h= g*u[i][i+1] - for (j=l; j < n; j++) - v[j][i]=u[i][j]/h - for (j=l; j < n; j++) - { - s=0.0 - for (k=l; k < n; k++) - s += u[i][k]*v[k][j] - for (k=l; k < n; k++) - v[k][j]+=(s*v[k][i]) - } - } - for (j=l; j < n; j++) - { - v[i][j] = 0; - v[j][i] = 0; - } - v[i][i] = 1; - g= e[i] - l= i - } - - // accumulation of left hand transformations - for (i=n-1; i != -1; i+= -1) - { - l= i+1 - g= q[i] - for (j=l; j < n; j++) - u[i][j] = 0; - if (g != 0.0) - { - h= u[i][i]*g - for (j=l; j < n; j++) - { - s=0.0 - for (k=l; k < m; k++) s += u[k][i]*u[k][j]; - f= s/h - for (k=i; k < m; k++) u[k][j]+=f*u[k][i]; - } - for (j=i; j < m; j++) u[j][i] = u[j][i]/g; - } - else - for (j=i; j < m; j++) u[j][i] = 0; - u[i][i] += 1; - } - - // diagonalization of the bidiagonal form - prec= prec*x - for (k=n-1; k != -1; k+= -1) - { - for (var iteration=0; iteration < itmax; iteration++) - { // test f splitting - var test_convergence = false - for (l=k; l != -1; l+= -1) - { - if (Math.abs(e[l]) <= prec) - { test_convergence= true - break - } - if (Math.abs(q[l-1]) <= prec) - break - } - if (!test_convergence) - { // cancellation of e[l] if l>0 - c= 0.0 - s= 1.0 - var l1= l-1 - for (i =l; i= itmax-1) - throw 'Error: no convergence.' - // shift from bottom 2x2 minor - x= q[l] - y= q[k-1] - g= e[k-1] - h= e[k] - f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y) - g= pythag(f,1.0) - if (f < 0.0) - f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x - else - f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x - // next QR transformation - c= 1.0 - s= 1.0 - for (i=l+1; i< k+1; i++) - { - g= e[i] - y= q[i] - h= s*g - g= c*g - z= pythag(f,h) - e[i-1]= z - c= f/z - s= h/z - f= x*c+g*s - g= -x*s+g*c - h= y*s - y= y*c - for (j=0; j < n; j++) - { - x= v[j][i-1] - z= v[j][i] - v[j][i-1] = x*c+z*s - v[j][i] = -x*s+z*c - } - z= pythag(f,h) - q[i-1]= z - c= f/z - s= h/z - f= c*g+s*y - x= -s*g+c*y - for (j=0; j < m; j++) - { - y= u[j][i-1] - z= u[j][i] - u[j][i-1] = y*c+z*s - u[j][i] = -y*s+z*c - } - } - e[l]= 0.0 - e[k]= f - q[k]= x - } - } - - //vt= transpose(v) - //return (u,q,vt) - for (i=0;i= 0; j--) - { - if (q[j] < q[i]) - { - // writeln(i,'-',j) - c = q[j] - q[j] = q[i] - q[i] = c - for(k=0;k b) + return a*Math.sqrt(1.0+(b*b/a/a)); + else if (b == 0.0) + return a; + return b*Math.sqrt(1.0+(a*a/b/b)); + } + + //Householder's reduction to bidiagonal form + + var f= 0.0; + var g= 0.0; + var h= 0.0; + var x= 0.0; + var y= 0.0; + var z= 0.0; + var s= 0.0; + + for (i=0; i < n; i++) + { + e[i]= g; + s= 0.0; + l= i+1; + for (j=i; j < m; j++) + s += (u[j][i]*u[j][i]); + if (s <= tolerance) + g= 0.0; + else + { + f= u[i][i]; + g= Math.sqrt(s); + if (f >= 0.0) g= -g; + h= f*g-s; + u[i][i]=f-g; + for (j=l; j < n; j++) + { + s= 0.0; + for (k=i; k < m; k++) + s += u[k][i]*u[k][j]; + f= s/h; + for (k=i; k < m; k++) + u[k][j]+=f*u[k][i]; + } + } + q[i]= g; + s= 0.0; + for (j=l; j < n; j++) + s= s + u[i][j]*u[i][j]; + if (s <= tolerance) + g= 0.0; + else + { + f= u[i][i+1]; + g= Math.sqrt(s); + if (f >= 0.0) g= -g; + h= f*g - s; + u[i][i+1] = f-g; + for (j=l; j < n; j++) e[j]= u[i][j]/h; + for (j=l; j < m; j++) + { + s=0.0; + for (k=l; k < n; k++) + s += (u[j][k]*u[i][k]); + for (k=l; k < n; k++) + u[j][k]+=s*e[k]; + } + } + y= Math.abs(q[i])+Math.abs(e[i]); + if (y>x) + x=y; + } + + // accumulation of right hand gtransformations + for (i=n-1; i != -1; i+= -1) + { + if (g != 0.0) + { + h= g*u[i][i+1]; + for (j=l; j < n; j++) + v[j][i]=u[i][j]/h; + for (j=l; j < n; j++) + { + s=0.0; + for (k=l; k < n; k++) + s += u[i][k]*v[k][j]; + for (k=l; k < n; k++) + v[k][j]+=(s*v[k][i]); + } + } + for (j=l; j < n; j++) + { + v[i][j] = 0; + v[j][i] = 0; + } + v[i][i] = 1; + g= e[i]; + l= i; + } + + // accumulation of left hand transformations + for (i=n-1; i != -1; i+= -1) + { + l= i+1; + g= q[i]; + for (j=l; j < n; j++) + u[i][j] = 0; + if (g != 0.0) + { + h= u[i][i]*g; + for (j=l; j < n; j++) + { + s=0.0; + for (k=l; k < m; k++) s += u[k][i]*u[k][j]; + f= s/h; + for (k=i; k < m; k++) u[k][j]+=f*u[k][i]; + } + for (j=i; j < m; j++) u[j][i] = u[j][i]/g; + } + else + for (j=i; j < m; j++) u[j][i] = 0; + u[i][i] += 1; + } + + // diagonalization of the bidiagonal form + prec= prec*x; + for (k=n-1; k != -1; k+= -1) + { + for (var iteration=0; iteration < itmax; iteration++) + { // test f splitting + var test_convergence = false; + for (l=k; l != -1; l+= -1) + { + if (Math.abs(e[l]) <= prec) + { test_convergence= true; + break; + } + if (Math.abs(q[l-1]) <= prec) + break; + } + if (!test_convergence) + { // cancellation of e[l] if l>0 + c= 0.0; + s= 1.0; + var l1= l-1; + for (i =l; i= itmax-1) + throw 'Error: no convergence.'; + // shift from bottom 2x2 minor + x= q[l]; + y= q[k-1]; + g= e[k-1]; + h= e[k]; + f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); + g= pythag(f,1.0); + if (f < 0.0) + f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x; + else + f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x; + // next QR transformation + c= 1.0; + s= 1.0; + for (i=l+1; i< k+1; i++) + { + g= e[i]; + y= q[i]; + h= s*g; + g= c*g; + z= pythag(f,h); + e[i-1]= z; + c= f/z; + s= h/z; + f= x*c+g*s; + g= -x*s+g*c; + h= y*s; + y= y*c; + for (j=0; j < n; j++) + { + x= v[j][i-1]; + z= v[j][i]; + v[j][i-1] = x*c+z*s; + v[j][i] = -x*s+z*c; + } + z= pythag(f,h); + q[i-1]= z; + c= f/z; + s= h/z; + f= c*g+s*y; + x= -s*g+c*y; + for (j=0; j < m; j++) + { + y= u[j][i-1]; + z= u[j][i]; + u[j][i-1] = y*c+z*s; + u[j][i] = -y*s+z*c; + } + } + e[l]= 0.0; + e[k]= f; + q[k]= x; + } + } + + //vt= transpose(v) + //return (u,q,vt) + for (i=0;i= 0; j--) + { + if (q[j] < q[i]) + { + // writeln(i,'-',j) + c = q[j]; + q[j] = q[i]; + q[i] = c; + for(k=0;k