BICPL: matrix_svd.c Source File

00001 /* ---------------------------------------------------------------------------- 00002 @COPYRIGHT : 00003 Copyright 1993,1994,1995 David MacDonald, 00004 McConnell Brain Imaging Centre, 00005 Montreal Neurological Institute, McGill University. 00006 Permission to use, copy, modify, and distribute this 00007 software and its documentation for any purpose and without 00008 fee is hereby granted, provided that the above copyright 00009 notice appear in all copies. The author and McGill University 00010 make no representations about the suitability of this 00011 software for any purpose. It is provided "as is" without 00012 express or implied warranty. 00013 ---------------------------------------------------------------------------- */ 00014 00015 #include <volume_io/internal_volume_io.h> 00016 #include <bicpl/numerical.h> 00017 00018 #ifndef lint 00019 static char rcsid[] = "$Header: /software/source//libraries/bicpl/Numerical/matrix_svd.c,v 1.6 2000/02/06 15:30:40 stever Exp $"; 00020 #endif 00021 00022 #define MAX_ITERATIONS 30 00023 00024 #define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \ 00025 (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0)) 00026 00027 #define TAKE_SIGN( a, b ) ((b) >= 0.0 ? fabs(a) : -fabs(a)) 00028 00029 /* ----------------------------- MNI Header ----------------------------------- 00030 @NAME : singular_value_decomposition 00031 @INPUT : m 00032 n 00033 a - array of size [m] by [n] 00034 @OUTPUT : w 00035 v 00036 @RETURNS : TRUE if successful 00037 @DESCRIPTION: Performs singular value decomposition of a matrix, based on 00038 the numerical recipes algorithms book. 00039 @METHOD : 00040 @GLOBALS : 00041 @CALLS : 00042 @CREATED : Jul 4, 1995 David MacDonald 00043 @MODIFIED : 00044 ---------------------------------------------------------------------------- */ 00045 00046 public BOOLEAN singular_value_decomposition( 00047 int m, 00048 int n, 00049 Real **a, 00050 Real w[], 00051 Real **v ) 00052 { 00053 BOOLEAN flag; 00054 int i, iter, j, jj, k, l, nm; 00055 Real c, f, h, s, x, y, z; 00056 Real anorm, g, scale, t; 00057 Real *rv1; 00058 Real at, bt, ct; 00059 00060 if( m < n ) 00061 { 00062 print_error("singular_value_decomposition(): matrix must have zero rows added to it.\n"); 00063 return( FALSE ); 00064 } 00065 00066 ALLOC( rv1, n ); 00067 00068 anorm = 0.0; 00069 scale = 0.0; 00070 g = 0.0; 00071 00072 for_less( i, 0, n ) 00073 { 00074 l = i + 1; 00075 rv1[i] = scale * g; 00076 g = 0.0; 00077 s = 0.0; 00078 scale = 0.0; 00079 00080 if( i < m ) 00081 { 00082 for_less( k, i, m ) 00083 scale += fabs( a[k][i] ); 00084 00085 if( scale != 0.0 ) 00086 { 00087 for_less( k, i, m ) 00088 { 00089 a[k][i] /= scale; 00090 s += a[k][i] * a[k][i]; 00091 } 00092 00093 f = a[i][i]; 00094 g = -TAKE_SIGN( sqrt(s), f ); 00095 h = f * g - s; 00096 a[i][i] = f - g; 00097 00098 if( i != n-1 ) 00099 { 00100 for_less( j, l, n ) 00101 { 00102 s = 0.0; 00103 for_less( k, i, m ) 00104 s += a[k][i] * a[k][j]; 00105 f = s / h; 00106 for_less( k, i, m ) 00107 a[k][j] += f * a[k][i]; 00108 } 00109 } 00110 00111 for_less( k, i, m ) 00112 a[k][i] *= scale; 00113 } 00114 } 00115 00116 w[i] = scale * g; 00117 g = 0.0; 00118 s = 0.0; 00119 scale = 0.0; 00120 if( i < m && i != n-1 ) 00121 { 00122 for_less( k, l, n ) 00123 scale += fabs( a[i][k] ); 00124 00125 if( scale != 0.0 ) 00126 { 00127 for_less( k, l, n ) 00128 { 00129 a[i][k] /= scale; 00130 s += a[i][k] * a[i][k]; 00131 } 00132 00133 f = a[i][l]; 00134 g = -TAKE_SIGN( sqrt(s), f ); 00135 h = f * g - s; 00136 a[i][l] = f - g; 00137 for_less( k, l, n ) 00138 rv1[k] = a[i][k] / h; 00139 if( i != m-1 ) 00140 { 00141 for_less( j, l, m ) 00142 { 00143 s = 0.0; 00144 for_less( k, l, n ) 00145 s += a[j][k] * a[i][k]; 00146 00147 for_less( k, l, n ) 00148 a[j][k] += s * rv1[k]; 00149 } 00150 } 00151 00152 for_less( k, l, n ) 00153 a[i][k] *= scale; 00154 } 00155 } 00156 00157 t = fabs(w[i] ) + fabs(rv1[i]); 00158 00159 if( t > anorm ) 00160 anorm = t; 00161 } 00162 00163 for_down( i, n-1, 0 ) 00164 { 00165 if( i < n-1 ) 00166 { 00167 if( g != 0.0 ) 00168 { 00169 for_less( j, l, n ) 00170 v[j][i] = (a[i][j] / a[i][l]) / g; 00171 00172 for_less( j, l, n ) 00173 { 00174 s = 0.0; 00175 for_less( k, l, n ) 00176 s += a[i][k] * v[k][j]; 00177 00178 for_less( k, l, n ) 00179 v[k][j] += s * v[k][i]; 00180 } 00181 } 00182 00183 for_less( j, l, n ) 00184 { 00185 v[i][j] = 0.0; 00186 v[j][i] = 0.0; 00187 } 00188 } 00189 00190 v[i][i] = 1.0; 00191 g = rv1[i]; 00192 l = i; 00193 } 00194 00195 for_down( i, n-1, 0 ) 00196 { 00197 l = i + 1; 00198 g = w[i]; 00199 if( i < n-1 ) 00200 { 00201 for_less( j, l, n ) 00202 a[i][j] = 0.0; 00203 } 00204 00205 if( g != 0.0 ) 00206 { 00207 g = 1.0 / g; 00208 if( i != n-1 ) 00209 { 00210 for_less( j, l, n ) 00211 { 00212 s = 0.0; 00213 for_less( k, l, m ) 00214 s += a[k][i] * a[k][j]; 00215 f = (s / a[i][i]) * g; 00216 for_less( k, i, m ) 00217 a[k][j] += f * a[k][i]; 00218 } 00219 } 00220 00221 for_less( j, i, m ) 00222 a[j][i] *= g; 00223 } 00224 else 00225 { 00226 for_less( j, i, m ) 00227 a[j][i] = 0.0; 00228 } 00229 00230 a[i][i] += 1.0; 00231 } 00232 00233 for_down( k, n-1, 0 ) 00234 { 00235 for_less( iter, 0, MAX_ITERATIONS ) 00236 { 00237 flag = TRUE; 00238 for_down( l, k, 0 ) 00239 { 00240 nm = l - 1; 00241 if( fabs(rv1[l]) + anorm == anorm ) 00242 { 00243 flag = FALSE; 00244 break; 00245 } 00246 00247 if( fabs(w[nm]) + anorm == anorm) 00248 break; 00249 } 00250 00251 if( flag ) 00252 { 00253 c = 0.0; 00254 s = 1.0; 00255 00256 for_less( i, l, k+1 ) 00257 { 00258 f = s * rv1[i]; 00259 00260 if( fabs(f)+anorm != anorm ) 00261 { 00262 g = w[i]; 00263 h = PYTHAG( f, g ); 00264 w[i] = h; 00265 h = 1.0 / h; 00266 c = g * h; 00267 s = -f * h; 00268 for_less( j, 0, m ) 00269 { 00270 y = a[j][nm]; 00271 z = a[j][i]; 00272 a[j][nm] = y*c + z*s; 00273 a[j][i] = z*c - y*s; 00274 } 00275 } 00276 } 00277 } 00278 00279 z = w[k]; 00280 if( l == k ) 00281 { 00282 if( z < 0.0 ) 00283 { 00284 w[k] = -z; 00285 for_less( j, 0, n ) 00286 v[j][k] = -v[j][k]; 00287 } 00288 break; 00289 } 00290 00291 if( iter == MAX_ITERATIONS-1 ) 00292 { 00293 print_error( "No convergence in %d SVD iterations.\n", 00294 MAX_ITERATIONS ); 00295 FREE( rv1 ); 00296 return( FALSE ); 00297 } 00298 00299 x = w[l]; 00300 nm = k - 1; 00301 y = w[nm]; 00302 g = rv1[nm]; 00303 h = rv1[k]; 00304 f = ((y-z)*(y+z)+(g-h)*(g+h)) / (2.0*h*y); 00305 g = PYTHAG( f, 1.0 ); 00306 f = ((x-z)*(x+z)+h*((y/(f+TAKE_SIGN(g,f)))-h)) / x; 00307 c =1.0; 00308 s =1.0; 00309 00310 for_less( j, l, nm+1 ) 00311 { 00312 i = j + 1; 00313 g = rv1[i]; 00314 y = w[i]; 00315 h = s * g; 00316 g = c * g; 00317 z = PYTHAG( f, h ); 00318 rv1[j] = z; 00319 c = f / z; 00320 s = h / z; 00321 f = x*c + g*s; 00322 g = g*c - x*s; 00323 h = y * s; 00324 y = y * c; 00325 for_less( jj, 0, n ) 00326 { 00327 x = v[jj][j]; 00328 z = v[jj][i]; 00329 v[jj][j] = x*c + z*s; 00330 v[jj][i] = z*c - x*s; 00331 } 00332 00333 z = PYTHAG( f, h ); 00334 w[j] = z; 00335 00336 if( z ) 00337 { 00338 z = 1.0 / z; 00339 c = f * z; 00340 s = h * z; 00341 } 00342 00343 f = c*g + s*y; 00344 x = c*y - s*g; 00345 00346 for_less( jj, 0, m ) 00347 { 00348 y = a[jj][j]; 00349 z = a[jj][i]; 00350 a[jj][j] = y*c + z*s; 00351 a[jj][i] = z*c - y*s; 00352 } 00353 } 00354 00355 rv1[l] = 0.0; 00356 rv1[k] = f; 00357 w[k] = x; 00358 } 00359 } 00360 00361 FREE( rv1 ); 00362 00363 return( TRUE ); 00364 }