BCLS: Bound Constrained Least Squares

Version 0.1

bclsqr.c

Go to the documentation of this file.
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