1 /* sv2uv.c CCMATH mathematics library source code.
2 *
3 * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
4 * This code may be redistributed under the terms of the GNU library
5 * public license (LGPL). ( See the lgpl.license file for details.)
6 * ------------------------------------------------------------------------
7 */
8 #include <stdlib.h>
9 #include <math.h>
10 #include "matutl.h"
sv2uv(double * d,double * a,double * u,int m,double * v,int n)11 int sv2uv(double *d,double *a,double *u,int m,double *v,int n)
12 { double *p,*p1,*q,*pp,*w,*e;
13 double s,t,h,r,sv;
14 int i,j,k,mm,nm,ms;
15 if(m<n) return -1;
16 w=(double *)calloc(m+n,sizeof(double)); e=w+m;
17 for(i=0,mm=m,p=a; i<n ;++i,--mm,p+=n+1){
18 if(mm>1){ sv=h=0.;
19 for(j=0,q=p,s=0.; j<mm ;++j,q+=n){
20 w[j]= *q; s+= *q* *q;
21 }
22 if(s>0.){
23 h=sqrt(s); if(*p<0.) h= -h;
24 s+= *p*h; s=1./s; t=1./(w[0]+=h);
25 sv=1.+fabs(*p/h);
26 for(k=1,ms=n-i; k<ms ;++k){
27 for(j=0,q=p+k,r=0.; j<mm ;q+=n) r+=w[j++]* *q;
28 r=r*s;
29 for(j=0,q=p+k; j<mm ;q+=n) *q-=r*w[j++];
30 }
31 for(j=1,q=p; j<mm ;) *(q+=n)=w[j++]*t;
32 }
33 *p=sv; d[i]= -h;
34 }
35 if(mm==1) d[i]= *p;
36 }
37 ldumat(a,u,m,n);
38 for(i=0,q=a; i<n ;++i){
39 for(j=0; j<n ;++j,++q){
40 if(j<i) *q=0.;
41 else if(j==i) *q=d[i];
42 }
43 }
44 for(i=0,mm=n,nm=n-1,p=a; i<n ;++i,--mm,--nm,p+=n+1){
45 if(i && mm>1){ sv=h=0.;
46 for(j=0,q=p,s=0.; j<mm ;++j,q+=n){
47 w[j]= *q; s+= *q* *q;
48 }
49 if(s>0.){
50 h=sqrt(s); if(*p<0.) h= -h;
51 s+= *p*h; s=1./s; t=1./(w[0]+=h);
52 sv=1.+fabs(*p/h);
53 for(k=1,ms=n-i; k<ms ;++k){
54 for(j=0,q=p+k,r=0.; j<mm ;q+=n) r+=w[j++]* *q;
55 for(j=0,q=p+k,r*=s; j<mm ;q+=n) *q-=r*w[j++];
56 }
57 for(k=0,p1=u+i; k<m ;++k,p1+=m){
58 for(j=0,q=p1,r=0.; j<mm ;) r+=w[j++]* *q++;
59 for(j=0,q=p1,r*=s; j<mm ;) *q++ -=r*w[j++];
60 }
61 }
62 *p=sv; d[i]= -h;
63 }
64 if(mm==1) d[i]= *p;
65 p1=p+1;
66 if(nm>1){ sv=h=0.;
67 for(j=0,q=p1,s=0.; j<nm ;++j,++q) s+= *q* *q;
68 if(s>0.){
69 h=sqrt(s); if(*p1<0.) h= -h;
70 sv=1.+fabs(*p1/h);
71 s+= *p1*h; s=1./s; t=1./(*p1+=h);
72 for(k=n,ms=n*(n-i); k<ms ;k+=n){
73 for(j=0,q=p1,pp=p1+k,r=0.; j<nm ;++j) r+= *q++ * *pp++;
74 for(j=0,q=p1,pp=p1+k,r*=s; j<nm ;++j) *pp++ -=r* *q++;
75 }
76 for(j=1,q=p1+1; j<nm ;++j) *q++ *=t;
77 }
78 *p1=sv; e[i]= -h;
79 }
80 if(nm==1) e[i]= *p1;
81 }
82 ldvmat(a,v,n);
83 qrbdv(d,e,u,m,v,n);
84 for(i=0; i<n ;++i){
85 if(d[i]<0.){ d[i]= -d[i];
86 for(j=0,p=v+i; j<n ;++j,p+=n) *p= - *p;
87 }
88 }
89 free(w);
90 return 0;
91 }
92