#
# ipscal package
# functions to minimize strain
# functions for low rank approximation
#
# version 0.0.1  only metric, one matrix
#

svdRect<-function(
	data,
	w=matrix(1,dim(data)[1],dim(data)[2]),
	p=2,
	init="rand",
	ithist=TRUE,
	itmax=100,
	eps=1e-6) {
n<-dim(data)[1]; m<-dim(data)[2]; itel<-1
nn<-1:n; mm<-1:m; pp<-1:p; wr<-rowSums(w); wc<-colSums(w)
if (is.list(init))
	{
	x<-init[[1]]; y<-init[[2]]
	}
else {
	x<-normalMat(n,p); y<-normalMat(m,p)
	}
c<-x%*%t(y); strain<-sum(w*(data-c)^2); pstrain=strain
repeat {
    for (i in nn) {
	    wi<-w[i,]; di<-data[i,]; c<-crossprod(y,wi*y)
		x[i,]<-solve(c,crossprod(y,wi*di))
		}
    for (j in mm) {
	    wj<-w[,j]; dj<-data[,j]; c<-crossprod(x,wj*x)
		y[j,]<-solve(c,crossprod(x,wj*dj))
		}
    c<-x%*%t(y); strain<-sum(w*(data-c)^2)		
	if (ithist) 	{
		cat("Iteration: ",formatC(itel,digits=6,width=6),
			"   SStress:   ",formatC(pstrain,digits=6,width=12,format="f"),
			" ==>",formatC(strain,digits=6,width=12,format="f"),"\n")
		}
	if (((pstrain-strain)<eps) || (itel==itmax)) break()
	pstrain<-strain; itel<-itel+1
}
return(list(x=x,y=y,c=c,strain=strain))
}