001 package numal.highprecision;
002
003 import numal.highprecision.*;
004
005 public class Linear_algebra extends Object {
006
007 /** Calculates the singular value decomposition UDV^T of a given matrix A.
008 * The matrix is first transformed to bidiagonal form by calling {@link #hshreabid}, the two
009 * transforming matrices are calculated by calling {@link #psttfmmat} and {@link #pretfmmat},
010 * and finally the singular value decomposition is calculated by {@link #qrisngvaldecbid}.
011 * Dependencies:
012 * {@link #hshreabid},
013 * {@link #psttfmmat},
014 * {@link #pretfmmat},
015 * {@link #qrisngvaldecbid}.
016 * <ul>
017 * <li><b>input:</b></li>
018 * <li>em[0] : the machine precision</li>
019 * <li>em[2] : the relative precision in the singular values</li>
020 * <li>em[4] : the maximal number of iterations to be performed</li>
021 * <li>em[6] : the minimal nonneglectable singular value</li>
022 * <li><b>output:</b></li>
023 * <li>em[1] : the infinity norm of the matrix</li>
024 * <li>em[3] : the maximal neglected superdiagonal element</li>
025 * <li>em[5] : the number of iterations performed</li>
026 * <li>em[7] : the numerical rank of the matrix; i.e. the number of singular values greater than or equal to em[6]</li>
027 * </ul>
028 *
029 * @param a <b>input</b>: the given mxn matrix <b>output</b>: the matrix U in the singular value decomposition UDV^T
030 * @param m the number of rows of the matrix a
031 * @param n the number of columns of a, <b>precondition:</b> n ≤ m
032 * @param val <b>output</b>: the singular values 1 to n
033 * @param v <b>output</b>: the transpose of matrix V in the singular value decomposition
034 * @param em input and output arameters mentioned in the list above
035 * @return number of singular values not found, i.e., a number not equal to zero if the number of iterations exceeds em[4]
036 *
037 */
038 public static int qrisngvaldec(double a[][], int m, int n, double val[], double v[][], double em[]) {
039 int i;
040 double b[] = new double[n+1];
041
042 hshreabid(a,m,n,val,b,em);
043 psttfmmat(a,n,v,b);
044 pretfmmat(a,m,n,val);
045 i=qrisngvaldecbid(val,b,m,n,a,v,em);
046 return i;
047 }
048
049 /** Reduces an m x n matrix <i>a</i> to a bidiagonal form <i>B</i>.
050 * Method originally from class numal.LinearAlgebra - Lau (2004)
051 * section 3.12.1.A Other transformations : To bidiagonal form - real matrices p.212.
052 * Dependencies: tammat, mattam, elmcol, elmrow
053 *
054 * @param a <b>input:</b> m by n matrix <b>output:</b> data concerning the premultiplying and postmultiplying matrices
055 * @param m number of rows
056 * @param n number of columns
057 * @param d <b>output:</b> vector of the diagonal of the bidiagonal matrix <i>B</i>
058 * @param b <b>output:</b> vector of the super diagonal of the bidiagonal matrix <i>B</i>
059 * @param em <b>output:</b>
060 */
061 public static void hshreabid(double a[][], int m, int n, double d[], double b[], double em[]) {
062 int i,j,i1;
063 double norm,machtol,w,s,f,g,h;
064
065 norm=0.0;
066 for (i=1; i<=m; i++) {
067 w=0.0;
068 for (j=1; j<=n; j++) w += Math.abs(a[i][j]);
069 if (w > norm) norm=w;
070 }
071 machtol=em[0]*norm;
072 em[1]=norm;
073 for (i=1; i<=n; i++) {
074 i1=i+1;
075 s=Basic.tammat(i1,m,i,i,a,a);
076 if (s < machtol)
077 d[i]=a[i][i];
078 else {
079 f=a[i][i];
080 s += f*f;
081 d[i] = g = (f < 0.0) ? Math.sqrt(s) : -Math.sqrt(s);
082 h=f*g-s;
083 a[i][i]=f-g;
084 for (j=i1; j<=n; j++)
085 Basic.elmcol(i,m,j,i,a,a,Basic.tammat(i,m,i,j,a,a)/h);
086 }
087 if (i < n) {
088 s=Basic.mattam(i1+1,n,i,i,a,a);
089 if (s < machtol)
090 b[i]=a[i][i1];
091 else {
092 f=a[i][i1];
093 s += f*f;
094 b[i] = g = (f < 0.0) ? Math.sqrt(s) : -Math.sqrt(s);
095 h=f*g-s;
096 a[i][i1]=f-g;
097 for (j=i1; j<=m; j++)
098 Basic.elmrow(i1,n,j,i,a,a,Basic.mattam(i1,n,i,j,a,a)/h);
099 }
100 }
101 }
102 }
103
104 public static void psttfmmat(double a[][], int n, double v[][], double b[]) {
105 int i,i1,j;
106 double h;
107
108 i1=n;
109 v[n][n]=1.0;
110 for (i=n-1; i>=1; i--) {
111 h=b[i]*a[i][i1];
112 if (h < 0.0) {
113 for (j=i1; j<=n; j++) v[j][i]=a[i][j]/h;
114 for (j=i1; j<=n; j++)
115 Basic.elmcol(i1,n,j,i,v,v,Basic.matmat(i1,n,i,j,a,v));
116 }
117 for (j=i1; j<=n; j++) v[i][j]=v[j][i]=0.0;
118 v[i][i]=1.0;
119 i1=i;
120 }
121 }
122
123 public static void pretfmmat(double a[][], int m, int n, double d[]) {
124 int i,i1,j;
125 double g,h;
126
127 for (i=n; i>=1; i--) {
128 i1=i+1;
129 g=d[i];
130 h=g*a[i][i];
131 for (j=i1; j<=n; j++) a[i][j]=0.0;
132 if (h < 0.0) {
133 for (j=i1; j<=n; j++)
134 Basic.elmcol(i,m,j,i,a,a,Basic.tammat(i1,m,i,j,a,a)/h);
135 for (j=i; j<=m; j++) a[j][i] /= g;
136 } else
137 for (j=i; j<=m; j++) a[j][i]=0.0;
138 a[i][i] += 1.0;
139 }
140 }
141
142 public static int qrisngvaldecbid(double d[], double b[], int m, int n, double u[][], double v[][], double em[]) {
143 int n0,n1,k,k1,i,i1,count,max,rnk;
144 double tol,bmax,z,x,y,g,h,f,c,s,min;
145
146 tol=em[2]*em[1];
147 count=0;
148 bmax=0.0;
149 max=(int) em[4];
150 min=em[6];
151 rnk=n0=n;
152 do {
153 k=n;
154 n1=n-1;
155 while (true) {
156 k--;
157 if (k <= 0) break;
158 if (Math.abs(b[k]) >= tol) {
159 if (Math.abs(d[k]) < tol) {
160 c=0.0;
161 s=1.0;
162 for (i=k; i<=n1; i++) {
163 f=s*b[i];
164 b[i] *= c;
165 i1=i+1;
166 if (Math.abs(f) < tol) break;
167 g=d[i1];
168 d[i1]=h=Math.sqrt(f*f+g*g);
169 c=g/h;
170 s = -f/h;
171 Basic.rotcol(1,m,k,i1,u,c,s);
172 }
173 break;
174 }
175 } else {
176 if (Math.abs(b[k]) > bmax) bmax=Math.abs(b[k]);
177 break;
178 }
179 }
180 if (k == n1) {
181 if (d[n] < 0.0) {
182 d[n] = -d[n];
183 for (i=1; i<=n0; i++) v[i][n] = -v[i][n];
184 }
185 if (d[n] <= min) rnk--;
186 n=n1;
187 } else {
188 count++;
189 if (count > max) break;
190 k1=k+1;
191 z=d[n];
192 x=d[k1];
193 y=d[n1];
194 g = (n1 == 1) ? 0.0 : b[n1-1];
195 h=b[n1];
196 f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
197 g=Math.sqrt(f*f+1.0);
198 f=((x-z)*(x+z)+h*(y/((f < 0.0) ? f-g : f+g)-h))/x;
199 c=s=1.0;
200 for (i=k1+1; i<=n; i++) {
201 i1=i-1;
202 g=b[i1];
203 y=d[i];
204 h=s*g;
205 g *= c;
206 z=Math.sqrt(f*f+h*h);
207 c=f/z;
208 s=h/z;
209 if (i1 != k1) b[i1-1]=z;
210 f=x*c+g*s;
211 g=g*c-x*s;
212 h=y*s;
213 y *= c;
214 Basic.rotcol(1,n0,i1,i,v,c,s);
215 d[i1]=z=Math.sqrt(f*f+h*h);
216 c=f/z;
217 s=h/z;
218 f=c*g+s*y;
219 x=c*y-s*g;
220 Basic.rotcol(1,m,i1,i,u,c,s);
221 }
222 b[n1]=f;
223 d[n]=x;
224 }
225 } while (n > 0);
226 em[3]=bmax;
227 em[5]=count;
228 em[7]=rnk;
229 return n;
230 }
231
232 }