sparseqr/sparseqr.hh

Go to the documentation of this file.

/*----------------------------------------------------------------------------
Sparse QR solver for linear systems of equations or least squares problems
Copyright (C) 1996, 1997, 1998 Thomas H. Robey

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Library General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Library General Public License for more details.

You should have received a copy of the GNU Library General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

Inquiries can be directed to trobey@arc.unm.edu or

Thomas H. Robey
925 Madison NE
Albuquerque, NM  87110
USA
----------------------------------------------------------------------------*/

#ifndef sparseqr_hh
#define sparseqr_hh

#include <iostream>
#include <cstddef>
#include <cassert>
#include <cfloat>

#ifdef DEBUG
#  define QRSTAT
#endif

namespace sparseqr {

  class Packed {
  public:
    int *col;
    double *value;
    unsigned size;
    Packed(unsigned cnt=0);
    void init(unsigned cnt=0);
    ~Packed();
    void Sort();
  };

  std::ostream& operator<<(std::ostream& os, const Packed& p);

  class J {
  public:
    int row1,row2;
    double c,s;
    J *next,*previous;
    void *operator new(size_t,double,double,int,int,J *);
  };

  std::ostream& operator<<(std::ostream& os, const J& j);

  class Qmatrix {
  public:
    J *qbegin,*qend;
    double *ApplyQ(double *);
    Qmatrix();
    ~Qmatrix();
  };

  std::ostream& operator<<(std::ostream& os, const Qmatrix& q);

  class Smatrix {
  public:
    int row,col;
    double value;
    Smatrix *right,*down;
    int match(Smatrix *);
    Smatrix *MoveDown(Smatrix **);
    Smatrix *MoveRight(Smatrix **);
    Smatrix *MarkEntry(Smatrix **sptr);
    void *operator new(size_t size,int r,int c,double v,Smatrix *rptr=NULL,
           Smatrix *dptr=NULL);
  };

  class QR {
  public:
    int nr,nc,*prow,*pcol;
    bool del_prow, del_pcol;
    unsigned *nrow;
#ifdef QRSTAT
    long int fillin,flops,ngivens;
#endif
    double ztol,rtol,*b;
    Packed *pmatrix;
    Smatrix **srow,**scol;
    Qmatrix *q;
    void BackSolve(int rank,double * x);
    int givens(int,int,int);
    QR(Packed *qrpmatrix,int qrnr,int qrnc,double *qrb=NULL,Qmatrix *qrq=NULL,
       int *qrprow=NULL,int *qrpcol=NULL,double qrztol=DBL_EPSILON,
       double qrrtol=0.0);
    ~QR();
    void RemoveEntry(int,int);
    Smatrix *RemoveEntry(int,int,Smatrix *,Smatrix *);
    int RemoveMarkedEntries(int);
    int sparseqr();
  };

} // namespace sparseqr

#endif // sparseqr_hh