################ (2021) Dr.W.G. Lindner, Leichlingen DE ### gsBox Gram-Schmidt box ################ for ortho/gonal/normalization proj(v,u) = dot(u,v)/dot(u,u)*u -- project v on u normalize(u) = u / abs(u) -- round to n decimals after point rnd(u,n) = do( test(number(n),Fixnum=n,Fixnum=4), float(floor(u*10^Fixnum+0.5)/10^Fixnum)) -- rnd(123.567898765543,3) GramSchmidt2(B)= do( O=zero(2,2), O[1]=B[1], O[2]=B[2] - proj(B[2],O[1]), O) -- result as rows GramSchmidt3(B)= do( O=zero(3,3), O[1]=B[1], O[2]=B[2] - proj(B[2],O[1]), O[3]=B[3] - proj(B[3],O[2]) - proj(B[3],O[1]), O) GramSchmidt(B) = do( Odim1 = dim(B,1), Odim2 = dim(B,2), O = zero(Odim1,Odim2), O[1]= B[1], for( k,2,Odim1, O[k] = B[k] - sum(j,1,k-1, proj(B[k],O[j]) ) ), O) -- OrthoNormalBasis ONB(B) = do(onb = B, odim = dim(B,1), for(k,1,odim, onb[k] = B[k] / abs(B[k])), onb) ## adapted from G. Weigt QR() = for(k,1,100, -- QR iteration to compute eigenvalues and EV's AT = transpose(A), ATo = GramSchmidt(AT), -- make AT orthogonal ATn = ONB(ATo), -- make ATo orthonormal U = transpose(ATn), A = dot(transpose(U),A,U), V = dot(V,U) )