test_linear.cc

COOOL：CWP面向对象最优化库（CWP Object Oriented Optimization Library） COOOL是C++类的一个集合

#include <iostream.h>
#include <stdlib.h>
#include <C++.hh>
#include <ObjFcns.hh>
#include <Optim.hh>

static int verbose        = 0;
static int modelsize      = 3;
static int itMax 	  = 100;
static double tol	  = 0.0000001;
static double upper 	=  1.;
static double lower 	=  -1.;

static char *opt_string = "S:A:N:M:T:L:U:V";

int main(int argc, char *argv[])
{
	int c,i, seed=0;
	int total = 1;

        
        while (( c = getopt(argc, argv, opt_string)) != -1)
        switch(c) {
	  case 'N':
	    modelsize = atoi(optarg);
	    break;
	  case 'M':
	    itMax = atoi(optarg);
	    break;
	  case 'T':
	    tol = atof(optarg);
	    cerr << "tolerance is "<< tol <<endl;
	    break;
	  case 'U':
	    upper = atof(optarg);
	    break;
	  case 'L':
	    lower = atof(optarg);
	    break;
	  case 'V':
	    ++verbose;
	    break;
	  default:
	    cerr << argv[0] << ": bad argument \n";
	    exit(1);
	}

	if( argc == 1) {
	    cerr << "Usage:  " << argv[0] << endl;
	    cerr << " -N number of unknowns" << endl;
	    cerr << " [-U upper bound of model parameters ]" << endl;
	    cerr << " [-L lower bound of model parameters ]" << endl;
	    cerr << " [-M maximum number of iterations]" << endl;
	    cerr << " [-T tolerance of the residue]" << endl;
	    cerr << " [-V verbose ]" << endl;
	    exit(-1);
	}


	/*************************
	  The section is for constructing a SpaseMatrix 
	  ***********************************************/

 	// read in the number of non-zero parameters in the matrix
	int nz;
	cin >> nz;

	// read the row, column indices and values of the non-zero elements
	Vector<int> icol(nz), irow(nz);
	Vector<double>  elem(nz);
	for (i=0; i< nz; i++) {
	  cin >> irow[i] >> icol[i] >> elem[i]; 
	}

	// obtain the size of the sparse matrix
	modelsize = icol.max()+1;
	int ndata = irow.max()+1;

	//construct the sparse matrix
	SpaMatrix<double>* jacob = new SpaMatrix<double>(elem, irow, icol);


	/***********End of Constructing Sparse Matrix**************/

	/**************************
	  Constructing a Dense Matrix
	  *****************************************/
	//DensMatrix<double>* jacob = new DensMatrix<double>(ndata,modelsize);
	//for (int irow=0; irow<ndata; irow++)
	//for (int icol=0; icol < modelsize; icol++)
	//(*jacob)[irow][icol] = 1./(irow+icol+1);
	/***********End of Constructing Dense Matrix**************/

	
	// Read the Right hand side
	Vector<double>* rhs = new Vector<double>(ndata);
	cin >> (*rhs);
	
	
	//Constrcut optimization objects
	LinearForward* lop = new LinearForward(jacob);
	/****IF USING ART, UNCOMMENT THE FOLLOWING LINE**********/
	// QuadraticOptima* opt = new ART(modelsize, lop, rhs, itMax, tol, verbose);
	/****IF USING SIRT, UNCOMMENT THE FOLLOWING LINE**********/
	// QuadraticOptima* opt = new SIRT(modelsize, lop, rhs, itMax, tol, verbose);

	/****IF USING IRLS, UNCOMMENT THE FOLLOWING LINE**********/
	QuadraticOptima* opt = new IterativeReweightedLS(modelsize, lop, rhs, 
							 1, 2,itMax, tol, tol, verbose);
	

	/****IF USING CGLS, UNCOMMENT THE FOLLOWING LINE**********/
	// QuadraticOptima* opt = new LSConjugateGradient(modelsize, lop, rhs,
	// itMax, tol, verbose);

	// construct the model space with UPPER and LOWER bound
  	Vector<double> v(modelsize), upBound(modelsize), lowBound(modelsize);
	upBound = upper;
	lowBound = lower;
	v = (double)0.0;
	
	// constructing the initial model without UPPER and LOWER bound
	Model<double> m(v);
	// constructing the initial model with UPPER and LOWER bound
	// Model<double> m(upBound,lowBound,v);


	// DO the optimization, output is in mOpt
	Model<double> mOpt(m);
	mOpt	=	opt->optimizer(m);
	cerr <<"The number of iterations: "<< opt->numIterations() << endl;

	// Output Results and relative information
	ofstream outfp("obtain.dat",ios::out);
	outfp << "The best model \n"<< mOpt <<endl;
	outfp << endl;
	outfp << endl;
	outfp <<"The number of iterations: "<< opt->numIterations() << endl;
	outfp <<"The first residue: " << opt->firstResidue() << endl;
	outfp <<"The final residue: " << opt->finalResidue() << endl;
	outfp <<"The residue history: "<<opt->allResidue() << endl;
}