BCLS: Bound Constrained Least Squares

Version 0.1

bccgls.c

Go to the documentation of this file.
00001 /* bccgls.c
00002    $Revision: 281 $ $Date: 2006-12-13 09:03:53 -0800 (Wed, 13 Dec 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 
00027 */
00034 #ifdef __APPLE__
00035   #include <vecLib/vecLib.h>
00036 #else
00037   #include <cblas.h>
00038 #endif
00039 #include <string.h>
00040 #include <assert.h>
00041 #include <float.h>
00042 
00043 #include "bcls.h"
00044 #include "bclib.h"
00045 #include "bccgls.h"
00046 
00066 static void
00067 aprod_free( BCLS *ls, const int mode, const int m, const int nFree,
00068             double dxFree[], double y[] )
00069 {
00070     int    *ix = ls->ix;
00071     double *dx = ls->dx;  // Used as workspace.
00072 
00073     assert( mode == BCLS_PROD_A || mode == BCLS_PROD_At );
00074 
00075     if (mode == BCLS_PROD_A) {
00076 
00077         bcls_scatter( nFree, ix, dxFree, dx );          // dx(ix) = dxFree.
00078         bcls_aprod( ls, BCLS_PROD_A, nFree, ix, dx, y );// y = A dx.
00079 
00080     } else {
00081         
00082         bcls_aprod( ls, BCLS_PROD_At, nFree, ix, dx, y );// dx = A' y.
00083         bcls_gather( nFree, ix, dx, dxFree );            // dxFree = dx(ix);
00084     }
00085 
00086     return;
00087 }
00088 
00128 static int
00129 bcls_cgls( BCLS *ls, int m, int n, int kmax, double tol, double damp,
00130            double c[], double x[], double r[], double s[], double p[],
00131            int *itns, double *opt, FILE *nout )
00132 {    
00133     int
00134         k = 0, info = 0, converged = 0, unstable = 0;
00135     double
00136         xmax = 0.0, resNE, Arnorm, Arnorm0,
00137         normp, norms, normx = 0.0, alpha, beta,
00138         gamma, gamma1, delta;
00139     const int
00140         damped = damp > 0.0;
00141     const double
00142         eps = DBL_EPSILON;
00143 
00144     // Initialize.
00145     bcls_dload( n, 0.0, x, 1 );                  // x  = 0
00146     aprod_free( ls, BCLS_PROD_At, m, n, s, r );  // s  = A'r
00147     if (c) 
00148         cblas_daxpy( n, 1.0, c, 1, s, 1 );       // s += c
00149     cblas_dcopy( n, s, 1, p, 1 );                // p  = s
00150     Arnorm0    = cblas_dnrm2( n, s, 1 );
00151     gamma      = Arnorm0 * Arnorm0;
00152 
00153     if (nout) {
00154         fprintf(nout,"\n CGLS:\n");
00155         fprintf(nout,"  Rows............. %8d\t", m);
00156         fprintf(nout,"  Columns.......... %8d\n", n);
00157         fprintf(nout,"  Optimality tol... %8.1e\t", tol);
00158         fprintf(nout,"  Iteration limit.. %8d\n", kmax);
00159         fprintf(nout,"  Damp............. %8.1e\t", damp);
00160         fprintf(nout,"  Linear term...... %8s\n",
00161                 (c) ? "yes" : "no" );
00162         fprintf(nout, "\n\n%5s %17s %17s %9s %10s\n",
00163                 "k", "x(1)","x(n)","normx","resNE");
00164         fprintf(nout, "%5d %17.10e %17.10e %9.2e %10.2e\n",
00165                 k, x[0], x[n-1], normx, 1.0 );
00166     }
00167     // -----------------------------------------------------------------
00168     // Main loop.
00169     // -----------------------------------------------------------------
00170     while ( k < kmax  &&  !converged  &&  !unstable ) {
00171 
00172         k++;
00173 
00174         aprod_free( ls, BCLS_PROD_A, m, n, p, s );   // s = A p
00175 
00176         norms = cblas_dnrm2( m, s, 1 );
00177         delta = norms * norms;
00178         if (damped) {
00179             normp  = cblas_dnrm2( n, p, 1 );
00180             delta += damp * normp * normp;
00181         }
00182         if (delta <= eps)
00183             delta =  eps;
00184 
00185         alpha = gamma / delta;
00186 
00187         cblas_daxpy( n,  alpha, p, 1, x, 1 );        // x  = x + alpha*p
00188         cblas_daxpy( m, -alpha, s, 1, r, 1 );        // r  = r - alpha*s
00189         aprod_free( ls, BCLS_PROD_At, m, n, s, r );  // s  = A'r
00190         if (damped)
00191             cblas_daxpy( n, -damp, x, 1, s, 1 );     // s += - damp*x
00192         if (c)
00193             cblas_daxpy( n,   1.0, c, 1, s, 1 );     // s += c
00194 
00195         Arnorm = cblas_dnrm2( n, s, 1 );
00196         gamma1 = gamma;
00197         gamma  = Arnorm * Arnorm;
00198         beta   = gamma / gamma1;
00199         
00200         cblas_dscal( n, beta, p, 1 );                // p  = beta*p
00201         cblas_daxpy( n, 1.0, s, 1, p, 1);            // p += s 
00202 
00203         // Check convergence.
00204         normx     = cblas_dnrm2( n, x, 1 );
00205         xmax      = fmax( xmax, normx );
00206         converged = Arnorm <= Arnorm0 * tol;
00207         unstable  = normx * tol >= 1.0;
00208 
00209         // Output.
00210         resNE = Arnorm / Arnorm0;
00211         if (nout)
00212             fprintf(nout, "%5d %17.10e %17.10e %9.2e %10.2e\n",
00213                     k, x[0], x[n-1], normx, resNE );
00214         
00215     } // while
00216 
00217     double shrink = normx / xmax;
00218     if (converged)
00219         info = 0;
00220     else if (k >= kmax)
00221         info = 1;
00222     else if (shrink <= sqrt(tol))
00223         info = 2;
00224 
00225     // Output.
00226     if (nout) fprintf(nout,"\n");
00227     *itns = k;
00228     *opt  = resNE;
00229     return info;
00230 }
00231 
00257 int
00258 bcls_newton_step_cgls( BCLS *ls, int m, int nFree, int ix[], double damp,
00259                        int itnLim, double tol, double dxFree[], double x[],
00260                        double c[], double r[], int *itns, double *opt )
00261 {
00262     int
00263         j;
00264     double
00265         *u;
00266     const double
00267         damp2 = damp*damp;
00268 
00269     if (c) {
00270         u = ls->wrk_u;
00271         bcls_gather( nFree, ix, c, u );       // u  =   c(ix).
00272         cblas_dscal( nFree, -1.0, u, 1 );     // u  = - c(ix).
00273         if ( damp > 0.0 )                     // u += - damp2 x(ix).
00274             for (j = 0; j < nFree; j++)
00275                 u[j] -= damp2 * x[ix[j]];
00276     }
00277     else
00278         u = NULL;
00279 
00280     return bcls_cgls( ls, m, nFree, itnLim, tol, damp2,
00281                       u, dxFree, r, ls->wrk_v, ls->wrk_w,
00282                       itns, opt, ls->minor_file );
00283 }

Generated on Sun Mar 4 22:50:03 2007 by Doxygen 1.5.1