pinv.hpp

Go to the documentation of this file.

/***************************************************************************
                        pinv.h -  description
                       -------------------------
    begin                : June 2006
    copyright            : (C) 2006 Erwin Aertbelien
    email                : firstname.lastname@mech.kuleuven.ac.be

 History (only major changes)( AUTHOR-Description ) :
 
 ***************************************************************************
 *   This library is free software; you can redistribute it and/or         *
 *   modify it under the terms of the GNU Lesser General Public            *
 *   License as published by the Free Software Foundation; either          *
 *   version 2.1 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     *
 *   Lesser General Public License for more details.                       *
 *                                                                         *
 *   You should have received a copy of the GNU Lesser General Public      *
 *   License along with this library; if not, write to the Free Software   *
 *   Foundation, Inc., 59 Temple Place,                                    *
 *   Suite 330, Boston, MA  02111-1307  USA                                *
 *                                                                         *
 ***************************************************************************/


#ifndef PINV_H
#define PINV_H

#include <math.h>
/*
inline void MatrixOutput(std::ostream& os,double* a,int stridea,int rows,int cols) {
        int i,j;
        for (i=0;i<rows;++i) {
                for (j=0;j<cols;++j) {
                        os << a[stridea*i+j];
                        if (j!=cols-1) os << "\t";
                }
                if (i!=rows-1) os << "\n";
        }
}
*/

class PseudoInverse {
        int maxsize;
        double * tmp;
    double* y2;
public:

    const double* A;

    int strideA;

    const double* Lx;

    const double* Ly;

    int m;

    int n;

    bool prepared;
        double * U;
        int strideU;
        double * V;
        int strideV;
        double * S;
    int *ndx;

public:
        PseudoInverse(int _maxsize);

    int smallestSV(double eps);

        /*int pinv(double* A,int strideA,int m,int n,double *Ai,int strideAi,
                int nrOfIts,double eps);
        **
         * Weighted PseudoInverse
         *   - t = A*q
         *   - The norm of t is given by t'*Wt*t
         *   - The norm of q is given by q'*Wq*q
         *   - Wt is POSITIVE SEMI-DEFINITE diagonal matrix
         *   - Wq is POSITIVE DEFINITE diagonal matrix
         *   - Lt is a diagonal matrix of size m x m, such that Wt = Lt'*Lt
         *   - Lq is a diagonal matrix of size n x n, such that Wq = Lq'*Lq
         *
         * \todo 
         *   check Lq for POS DEF. COND
         *
        void wpinv(
                double* A,int strideA,int m,int n,double *Awi,int strideAwi,
            double* Lt,double* Lq,int nrOfIts,double eps);
    */
        ~PseudoInverse();

    int prepare( const double* _A,int _strideA,int _m,int _n,
                                 const double *_Ly,const double* _Lx,
                                 int nrOfIts);


    void inverseWithNullSpace(
        const double* y,
        double* x,
        const double* xd,
        double eps);
    void inverse(
        const double* y,
        double* x,
        double eps);

    double derivsv(
        const double* Adot, 
        int strideAdot,
        int svnr);

};



#endif