1 /* sv2val.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"
sv2val(double * d,double * a,int m,int n)11 int sv2val(double *d,double *a,int m,int n)
12 { double *p,*p1,*q,*w,*v;
13 double s,h,u;
14 int i,j,k,mm,nm,ms;
15 if(m<n) return -1;
16 w=(double *)calloc(m,sizeof(double));
17 for(i=0,mm=m,p=a; i<n && mm>1 ;++i,--mm,p+=n+1){
18 for(j=0,q=p,s=0.; j<mm ;++j,q+=n){
19 w[j]= *q; s+= *q* *q;
20 }
21 if(s>0.){
22 h=sqrt(s); if(*p<0.) h= -h;
23 s+= *p*h; s=1./s; w[0]+=h;
24 for(k=1,ms=n-i; k<ms ;++k){
25 for(j=0,q=p+k,u=0.; j<mm ;q+=n) u+=w[j++]* *q;
26 u=u*s;
27 for(j=0,q=p+k; j<mm ;q+=n) *q-=u*w[j++];
28 }
29 *p= -h;
30 }
31 }
32 for(i=0,p=a; i<n ;++i,p+=n){
33 for(j=0,q=p; j<i ;++j) *q++ =0.;
34 }
35 for(i=0,mm=n,nm=n-1,p=a; i<n ;++i,--mm,--nm,p+=n+1){
36 if(i && mm>1){
37 for(j=0,q=p,s=0.; j<mm ;++j,q+=n){
38 w[j]= *q; s+= *q* *q;
39 }
40 if(s>0.){
41 h=sqrt(s); if(*p<0.) h= -h;
42 s+= *p*h; s=1./s; w[0]+=h;
43 for(k=1,ms=n-i; k<ms ;++k){
44 for(j=0,q=p+k,u=0.; j<mm ;q+=n) u+=w[j++]* *q;
45 u*=s;
46 for(j=0,q=p+k; j<mm ;q+=n) *q-=u*w[j++];
47 }
48 *p= -h;
49 }
50 }
51 p1=p+1;
52 if(nm>1){
53 for(j=0,q=p1,s=0.; j<nm ;++j,++q) s+= *q* *q;
54 if(s>0.){
55 h=sqrt(s); if(*p1<0.) h= -h;
56 s+= *p1*h; s=1./s; *p1+=h;
57 for(k=n,ms=n*(m-i); k<ms ;k+=n){
58 for(j=0,q=p1,v=p1+k,u=0.; j<nm ;++j) u+= *q++ * *v++;
59 u*=s;
60 for(j=0,q=p1,v=p1+k; j<nm ;++j) *v++ -=u* *q++;
61 }
62 *p1= -h;
63 }
64 }
65 }
66 for(j=0,p=a; j<n ;++j,p+=n+1){
67 d[j]= *p; if(j<n-1) w[j]= *(p+1); else w[j]=0.;
68 }
69 qrbdi(d,w,n);
70 for(i=0; i<n ;++i) if(d[i]<0.) d[i]= -d[i];
71 free(w);
72 return 0;
73 }
74