---

F4ThomasLS.hpp

#ifndef __F4ThomasLS_H__
#define __F4ThomasLS_H__



#include "AbsF4TriDiagonalLS.hpp"

#include <math.h>
#include "sseUtil.h"


class F4ThomasLS : public AbsF4TriDiagonalLS
{
private:
                                                
Float4Vector*                   uOut;   
Float4Vector*                   dOut;
Float4Vector*                   lOut;   

Float4Vector*                   y;

protected:

// for parallel executions
inline void computeLRdecomposition()
        {
/*      
        for(int line=0; line<4; line++)
                {
                dOut->set0(line, 0, A->getD0(line, 0) );

                for(int j=1; j<N; j++)
                        {
                        lOut->set(line, j, A->getL(line, j)/dOut->get(line, j) );
        
                        dOut->set0(line, j, A->getD0(line, j) - lOut->get(line, j)*A->getU(line,j) );
                        }
                }
*/
        float *dout, *ad0, *al, *au, *lout;
//      int j4=4;


        dout= dOut->get0();
        ad0= A->getD0();
        al= A->getL();
        au= A->getU();
        lout= lOut->get(); 


/*      prefetcht0(al[j4]);                             // use one time
        prefetcht0(ad0[j4]);                    // use one time
        prefetcht0(au[j4]);                             // use one time
*/
        movaps_m2r(ad0[0], xmm0);               // dout0
        movaps_r2m(xmm0, dout[0]);      

        for(int j=4; j<N4; j+=4) // , j4+=4)
                {
/*              prefetcht0(al[j4]);                             // use one time
                prefetcht0(ad0[j4]);                            // use one time
                prefetcht0(au[j4]);                             // use one time
*/
                movaps_m2r(al[j], xmm1);                // al
//              divps_r2r(xmm0, xmm1);                  // al/dout
                rcpps_r2r(xmm0, xmm3);                  // dout -> 1/dout  division optimisation
                mulps_r2r(xmm3, xmm1);                  // al*(1/dout)
                movaps_r2m(xmm1, lout[j]);              // lout

                movaps_m2r(ad0[j], xmm0);               // ad0
                movaps_m2r(au[j], xmm2);                // au
                mulps_r2r(xmm2, xmm1);                  // lout*au
                subps_r2r(xmm1, xmm0);                  // ad0 - lout*au
                movaps_r2m(xmm0, dout[j]);
                }
        }

inline void forwardSubstitution()
        {
/*      
        for(int line=0; line<4; line++)
                {
                y->set0(line, 0, b->get0(line, 0) );

                for(int j=1; j<N; j++)
                        y->set0(line, j, b->get0(line, j) - lOut->get(line, j)*y->get(line, j) );
                }
*/
        float *yv, *bv, *lout;
//      int j4=4;

        yv= y->get0();
        bv= b->get0();
        lout= lOut->get();

/*      prefetcht0(lout[j4]);                   // use one time
        prefetcht0(bv[j4]);                             // use one time
*/
        movaps_m2r(bv[0], xmm0);                // b0
        movaps_r2m(xmm0, yv[0]);                // y0

        for(int j=4; j<N4; j+=4)  // , j4+=4)
                {
/*              prefetcht0(lout[j4]);                   // use one time
                prefetcht0(bv[j4]);                     // use one time
*/
                movaps_m2r(bv[j], xmm1);        // bj
                movaps_m2r(lout[j], xmm2);      // lout
                mulps_r2r(xmm2, xmm0);          // lout*y
                subps_r2r(xmm0, xmm1);          // bj - lout*y
                movaps_r2r(xmm1, xmm0);
                movaps_r2m(xmm0, yv[j]);        // yj
                }
        }

inline void backwardSubstitution()
        {
/*      for(int line=0; line<4; line++)
                {
                x->set(line, N, y->get(line, N)/dOut->get(line, N) );

                for(int j=N-1; j>0; j--)
                        x->set(line, j, (y->get(line,j) - A->getU(line,j)*x->get0(line,j))/dOut->get(line,j) );
                }
*/
        float *xv, *yv, *dout, *au;
//      int j4=N4-4;

        xv= x->get();
        yv= y->get();
        dout= dOut->get();
        au= A->getU();

        prefetcht1(au[N4-4]);                   // use one time
/*      prefetcht0(au[j4]);                     // use one time
        prefetcht0(dout[j4]);                   // use one time
        prefetcht0(yv[j4]);                             // use one time
*/
        movaps_m2r(yv[N4], xmm0);               // y
        movaps_m2r(dout[N4], xmm1);             // dout
//      divps_r2r(xmm1, xmm0);                  // y/dout
        rcpps_r2r(xmm1, xmm4);                  // dout -> 1/dout  division optimisation
        mulps_r2r(xmm4, xmm0);                  // y*(1/dout)

        movaps_r2m(xmm0, xv[N4]);               // x

        for(int j=N4-4; j>0; j -=4) // , j4-=4)
                {
                prefetcht1(au[j-4]);                    // use one time
/*              prefetcht0(au[j4]);                     // use one time
                prefetcht0(dout[j4]);                   // use one time
                prefetcht0(yv[j4]);                     // use one time
*/
                movaps_m2r(yv[j], xmm1);        // y
                movaps_m2r(au[j], xmm2);        // au
                movaps_m2r(dout[j], xmm3);      // dout

                mulps_r2r(xmm0, xmm2);          // x*au
                subps_r2r(xmm2, xmm1);          // y - x*au
//              divps_r2r(xmm3, xmm1);          // (y - x*au)/dout
                rcpps_r2r(xmm3, xmm4);                  // dout -> 1/dout  division optimisation
                mulps_r2r(xmm4, xmm1);          // (y - x*au)*(1/dout)
        
                movaps_r2r(xmm1, xmm0);
                movaps_r2m(xmm0, xv[j]);        // x
                }
        }               

// combined LR/forward  reducing data flow
// for parallel executions
inline void computeLR_forward()
        {
/*      
        for(int line=0; line<4; line++)
                {
                dOut->set0(line, 0, A->getD0(line, 0) );
                y->set0(line, 0, b->get0(line, 0) );

                for(int j=1; j<N; j++)
                        {
                        lOut->set(line, j, A->getL(line, j)/dOut->get(line, j) );
                        dOut->set0(line, j, A->getD0(line, j) - lOut->get(line, j)*A->getU(line,j) );

                        y->set0(line, j, b->get0(line, j) - lOut->get(line, j)*y->get(line, j) );
                        }
                }

*/
        float *dout, *ad0, *al, *au,  *yv, *bv;  // lout is not needed
//      int j4=4;

        dout= dOut->get0();
        ad0= A->getD0();
        al= A->getL();
        au= A->getU();
//      lout= lOut->get();                              // not needed 
        yv= y->get0();
        bv= b->get0();

        prefetcht1(bv[4]);                              // use one time
        prefetcht1(ad0[4]);                             // use one time
/*      prefetcht0(al[j4]);                             // use one time
        prefetcht0(ad0[j4]);                    // use one time
        prefetcht0(au[j4]);                             // use one time
        prefetcht0(bv[j4]);                             // use one time
*/
        movaps_m2r(ad0[0], xmm0);               // dout0
        movaps_r2m(xmm0, dout[0]);      

        movaps_m2r(bv[0], xmm4);                // b0
        movaps_r2m(xmm4, yv[0]);                // y0

        for(int j=4; j<N4; j+=4) // , j4+=4)
                {
                prefetcht1(bv[j+4]);                            // use one time
                prefetcht1(ad0[j+4]);                           // use one time
                // load next data in Lx
/*              prefetcht0(al[j4]);                     // use one time
                prefetcht0(ad0[j4]);                    // use one time
                prefetcht0(au[j4]);                     // use one time
                prefetcht0(bv[j4]);                     // use one time
*/
                movaps_m2r(bv[j], xmm5);        // bj

                movaps_m2r(al[j], xmm1);                // al
//              divps_r2r(xmm0, xmm1);                  // al/dout
                rcpps_r2r(xmm0, xmm3);                  // dout -> 1/dout  division optimisation
                mulps_r2r(xmm3, xmm1);                  // al*(1/dout)
//              movaps_r2m(xmm1, lout[j]);              // lout                         NOT NEEDED anymore in results

                movaps_m2r(ad0[j], xmm0);               // ad0
                movaps_m2r(au[j], xmm2);                // au
                mulps_r2r(xmm1, xmm2);                  // lout*au
                subps_r2r(xmm2, xmm0);                  // ad0 - lout*au
                movaps_r2m(xmm0, dout[j]);

                mulps_r2r(xmm1, xmm4);          // lout*y
                subps_r2r(xmm4, xmm5);          // bj - lout*y
                movaps_r2r(xmm5, xmm4);
                movaps_r2m(xmm4, yv[j]);        // yj
                }
        }

public:
        
F4ThomasLS(int asize):AbsF4TriDiagonalLS(asize)
        {
        uOut= new Float4Vector(N-1);
        dOut= new Float4Vector(N);
        lOut= new Float4Vector(N-1);

        y= new Float4Vector(N);
        }

F4ThomasLS(int asize, int lines):AbsF4TriDiagonalLS(asize, lines)
        {
        uOut= new Float4Vector(N-1);
        dOut= new Float4Vector(N);
        lOut= new Float4Vector(N-1);

        y= new Float4Vector(N);
        }

virtual ~F4ThomasLS()
        {
        delete uOut; 
        delete dOut; 
        delete lOut; 
        delete y;
        }

virtual void loadA(F4TriDiagonalMatrix* m);
virtual void loadMatrix(Float4Vector* u1, Float4Vector* d1, Float4Vector* l1);
virtual void loadX(Float4Vector* data);
virtual void loadB(Float4Vector* data);


inline void solve()
        {
//      computeLRdecomposition();
//      forwardSubstitution();

        computeLR_forward();                    // combined optimisation
        backwardSubstitution();
        }
        
virtual void output()
        {
        printf("F4Thomas ");
        AbsF4TriDiagonalLS::output();
        }

virtual void output(FILE* file)
        {
        fprintf(file, "F4Thomas ");
        AbsF4TriDiagonalLS::output(file);
        }

};

#endif

         

Restoreinpaint sourceforge project `C++/Java Image Processing, Restoration, Inpainting Project'.

Bernard De Cuyper: Open Project Leader: Concept, design and development. Bernard De Cuyper & Eddy Fraiha 2002, 2003. Bernard De Cuyper 2004. Open and free, for friendly usage only. Modifications on Belgium ground of this piece of artistic work, by governement institutions or companies, must be notified to Bernard De Cuyper.