Version 0.1
00001 /* bclsqr.c 00002 $Revision: 273 $ $Date: 2006-09-04 15:59:04 -0700 (Mon, 04 Sep 2006) $ 00003 00004 ---------------------------------------------------------------------- 00005 This file is part of BCLS (Bound-Constrained Least Squares). 00006 00007 Copyright (C) 2006 Michael P. Friedlander, Department of Computer 00008 Science, University of British Columbia, Canada. All rights 00009 reserved. E-mail: <mpf@cs.ubc.ca>. 00010 00011 BCLS is free software; you can redistribute it and/or modify it 00012 under the terms of the GNU Lesser General Public License as 00013 published by the Free Software Foundation; either version 2.1 of the 00014 License, or (at your option) any later version. 00015 00016 BCLS is distributed in the hope that it will be useful, but WITHOUT 00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 00018 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General 00019 Public License for more details. 00020 00021 You should have received a copy of the GNU Lesser General Public 00022 License along with BCLS; if not, write to the Free Software 00023 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 00024 USA 00025 ---------------------------------------------------------------------- 00026 */ 00032 #ifdef __APPLE__ 00033 #include <vecLib/vecLib.h> 00034 #else 00035 #include <cblas.h> 00036 #endif 00037 #include <string.h> 00038 #include <stdio.h> 00039 00040 #include "bcls.h" 00041 #include "bclib.h" 00042 #include "bclsqr.h" 00043 #include "lsqr.h" 00044 00078 static void 00079 aprod_free_lsqr( const int mode, const int mSubProb, const int nFree, 00080 double dxFree[], double y[], void *UsrWrk) 00081 { 00082 int j; 00083 BCLS *ls = (BCLS *)UsrWrk; // Reclaim access to the BCLS workspace. 00084 const int m = ls->m; 00085 const int preconditioned = ls->Usolve != NULL; 00086 const int damped = ls->damp_actual > 0.0; 00087 const double damp = ls->damp_actual; 00088 int *ix = ls->ix; 00089 double *dx = ls->dx; // Used as workspace. 00090 double *dy = ls->wrk_u; // ... 00091 00092 if (mode == 1) { 00093 00094 // Solve U dy = dxFree. U is nFree-by-nFree, and so the 00095 // relevant part of dy has length nFree. 00096 if (preconditioned) 00097 bcls_usolve( ls, BCLS_PRECON_U, nFree, ix, dy, dxFree ); 00098 else 00099 cblas_dcopy( nFree, dxFree, 1, dy, 1 ); 00100 00101 // y2 <- y2 + damp * dy(1:nFree). 00102 if (damped) 00103 cblas_daxpy( nFree, damp, dy, 1, &y[m], 1 ); 00104 00105 // Scatter dy into dx: dx(ix) <- dy(1:nFree). 00106 for (j = 0; j < nFree; j++) dx[ ix[j] ] = dy[j]; 00107 00108 // dy <- A(:,ix) * dx(ix). 00109 bcls_aprod( ls, BCLS_PROD_A, nFree, ix, dx, dy ); 00110 00111 // y1 <- y1 + dy. 00112 cblas_daxpy( m, 1.0, dy, 1, y, 1 ); 00113 00114 } 00115 else { // mode == 2 00116 00117 // dx <- A' * y1. 00118 bcls_aprod( ls, BCLS_PROD_At, nFree, ix, dx, y ); 00119 00120 // Gather dx(ix) into dy: dy(1:nFree) <- dx(ix). 00121 for (j = 0; j < nFree; j++) dy[j] = dx[ ix[j] ]; 00122 00123 // dy <- dy + damp * y2. 00124 if (damped) 00125 cblas_daxpy( nFree, damp, &y[m], 1, dy, 1 ); 00126 00127 // Solve U' dx = dy. 00128 if (preconditioned) 00129 bcls_usolve( ls, BCLS_PRECON_Ut, nFree, ix, dy, dx ); 00130 else 00131 cblas_dcopy( nFree, dy, 1, dx, 1 ); 00132 00133 // dxFree <- dxFree + dx. 00134 cblas_daxpy( nFree, 1.0, dx, 1, dxFree, 1 ); 00135 00136 } 00137 return; 00138 } 00139 00167 int 00168 bcls_newton_step_lsqr( BCLS *ls, int m, int nFree, int ix[], double damp, 00169 int itnLim, double tol, double dxFree[], double x[], 00170 double c[], double r[], int *itns, double *opt ) 00171 { 00172 int j, k; // Misc. counters. 00173 int mpn; // No. of rows in [ N; damp I ]. 00174 int unscale_dxFree = 0; 00175 const int rescaling_method = 1; 00176 const int linear = c != NULL; 00177 const int damped = damp > 0.0; 00178 const double damp_min = ls->damp_min; 00179 const double damp2 = damp * damp; 00180 const double zero = 0.0; 00181 double damp_actual; 00182 double beta1, beta2; 00183 00184 // LSQR outputs 00185 int istop; // Termination flag. 00186 double anorm; // Estimate of Frobenious norm of Abar. 00187 double acond; // Estimate of condition no. of Abar. 00188 double rnorm; // Estimate of the final value of norm(rbar). 00189 double xnorm; // Estimate of the norm of the final solution dx. 00190 00191 // Set r(m+1:) <- - 1/beta c + damp^2/beta x, where 00192 // beta = max(min_damp, damp). 00193 // Note that r is declared length (m+n), so there is always enough 00194 // space. The chosen damping parameter is stored in 00195 // ls->damp_actual so that it can be used in bcls_aprod_free. 00196 if (!damped && !linear) { 00197 mpn = m; 00198 ls->damp_actual = 0.0; 00199 } 00200 else { 00201 00202 mpn = m + nFree; 00203 00204 //-------------------------------------------------------------- 00205 // Rescale the RHS. Will need to unscale LSQR's solution later. 00206 //-------------------------------------------------------------- 00207 if (rescaling_method) { 00208 00209 if (linear) { 00210 00211 unscale_dxFree = 1; 00212 damp_actual = fmax( damp, damp_min ); 00213 ls->damp_actual = damp_actual; 00214 00215 // r(1:m) = damp_actual * r(1:m) 00216 cblas_dscal( m, damp_actual, r, 1 ); 00217 00218 // r(m+1:) = - c 00219 for (j = 0; j < nFree; j++ ) 00220 r[m+j] = - c[ ix[j] ]; 00221 00222 // r(m+1:) = - c - damp^2 * x 00223 if (damped) 00224 for (j = 0; j < nFree; j++) 00225 r[m+j] -= damp2 * x[ ix[j] ]; 00226 } 00227 else { 00228 00229 ls->damp_actual = damp; 00230 for (j = 0; j < nFree; j++) 00231 r[m+j] = - damp * x[ ix[j] ]; 00232 } 00233 } 00234 //-------------------------------------------------------------- 00235 // No rescaling. 00236 //-------------------------------------------------------------- 00237 else { 00238 if (damped && linear) { 00239 00240 damp_actual = fmax( damp, damp_min ); 00241 beta1 = -1.0 / damp_actual; 00242 beta2 = -damp * damp / damp_actual; 00243 ls->damp_actual = damp_actual; 00244 00245 for (j = 0; j < nFree; j++) { 00246 k = ix[j]; 00247 r[m+j] = beta1 * c[k] + beta2 * x[k]; 00248 } 00249 } 00250 else if ( damped && !linear ) { 00251 00252 beta2 = - damp; 00253 ls->damp_actual = damp; 00254 for (j = 0; j < nFree; j++) 00255 r[m+j] = beta2 * x[ ix[j] ]; 00256 } 00257 else if (!damped && linear ) { 00258 00259 beta1 = - 1.0 / damp_min; 00260 ls->damp_actual = damp_min; 00261 for (j = 0; j < nFree; j++) 00262 r[m+j] = beta1 * c[ ix[j] ]; 00263 } 00264 } 00265 } 00266 00267 // ----------------------------------------------------------------- 00268 // Solve the subproblem with LSQR. 00269 // ----------------------------------------------------------------- 00270 bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR]), BCLS_TIMER_START ); 00271 00272 lsqr( mpn, nFree, aprod_free_lsqr, zero, (void *)ls, 00273 r, ls->wrk_v, ls->wrk_w, dxFree, NULL, 00274 tol, tol, ls->conlim, itnLim, ls->minor_file, 00275 &istop, itns, &anorm, &acond, 00276 &rnorm, opt, &xnorm ); 00277 00278 bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR]), BCLS_TIMER_STOP ); 00279 00280 // ----------------------------------------------------------------- 00281 // Cleanup and exit. First unscale dxFree if needed. 00282 // ----------------------------------------------------------------- 00283 if (rescaling_method && unscale_dxFree) 00284 cblas_dscal( nFree, 1.0/damp_actual, dxFree, 1 ); 00285 00286 if (istop <= 3) 00287 return 0; 00288 else if (istop == 5) 00289 return 1; 00290 else 00291 return 2; 00292 }
Generated on Sun Mar 4 22:50:03 2007 by Doxygen 1.5.1