/*
 $Id: block.cc,v 1.2 1996/10/04 15:06:24 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "block.h"
#include "utils.h"

#include "numerics.h"
#include "dirichlet.h"
#include "family.h"

//-------------------------------------------------------------------------


//	Multiplication routines for Blocks:

//	 LMult: multiply with lower Block
//	 UMult:		      upper Block

//	 LTMult:              transposed lower Block
//	 UTMult:              transposed upper Block

//-------------------------------------------------------------------------

StaticAllocator<BlockNode> 	    BlockNode:: 	 alloc(100);
StaticAllocator<NeighbourBlock>     NeighbourBlock:: 	 alloc(100);
StaticAllocator<AsymNeighbourBlock> AsymNeighbourBlock:: alloc(100);

//-------------------------------------------------------------------------


void BlockNode:: reset(int nComp)
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i)  D[i] = 0.0;
}
//-------------------------------------------------------------------------

void NeighbourBlock:: reset(int nComp)
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i)  L[i] = 0.0;
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: reset(int nComp)
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i)  L[i] = U[i] = 0.0;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void BlockNode:: setTranspose(int nComp)
{
    int i, k;
    Vector<Num*> Ap(nComp);

    Ap[1] = D;
    for (i=2; i<=nComp; ++i)  Ap[i] = Ap[i-1] + nComp;

    FORNCOMP(i) for (k=1; k<i; ++k)  Ap[k][i] = Ap[i][k];
}
//-------------------------------------------------------------------------


void BlockNode:: Invert(Matrix<Num>& AM, Matrix<Num>& AMInv)
{
    int i, j, n;
    const int nComp = AM.cHigh();

    n = 1;
    FORNCOMP(i) FORNCOMP(j) AM(i,j) = D[n++];

    ::invert(AM, AMInv);

    n = 1;
    FORNCOMP(i) FORNCOMP(j) DInv[n++] = AMInv(i,j);
}
//-------------------------------------------------------------------------

void BlockNode:: Add(Vector<Num*>& Ap)
{
    int i,k,n=1;
    const int nComp = Ap.high();
    FORNCOMP(i) FORNCOMP(k) D[n++] += Ap[i][k];
}
//-------------------------------------------------------------------------

void BlockNode:: getD(Vector<Num>& data, int nComp) const
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i) data[i] = D[i];
}
//-------------------------------------------------------------------------

void BlockNode:: Mult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[i] += D[n++] * rhs[k];
}
//-------------------------------------------------------------------------

void BlockNode:: TMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k) lhs[k] += conj(D[n++]) * rhs[i];
}
//-------------------------------------------------------------------------

void BlockNode:: Div(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k) lhs[i] += DInv[n++] * rhs[k];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void NeighbourBlock:: LAdd(Vector<Num*>& Ap)
{
    int i,k,n=1;
    const int nComp = Ap.high();
    FORNCOMP(i) FORNCOMP(k) L[n++] += Ap[i][k];
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: UAdd(Vector<Num*>& Ap)
{
    int i,k,n=1;
    const int nComp = Ap.high();
    FORNCOMP(i) FORNCOMP(k) U[n++] += Ap[i][k];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void NeighbourBlock:: getL(Vector<Num>& data, int nComp) const
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i) data[i] = L[i];
}
//-------------------------------------------------------------------------

void NeighbourBlock:: getU(Vector<Num>& data, int nComp) const
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i) data[i] = L[i];
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: getU(Vector<Num>& data, int nComp) const
{
    int i;
    const int n = nComp*nComp;
    for (i=1; i<=n; ++i) data[i] = U[i];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void NeighbourBlock:: LMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[i] += L[n++] * rhs[k];
}
//-------------------------------------------------------------------------

void NeighbourBlock:: UMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[k] += L[n++] * rhs[i];
}
//-------------------------------------------------------------------------

void NeighbourBlock:: LTMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[k] += conj(L[n++]) * rhs[i];
}
//-------------------------------------------------------------------------

void NeighbourBlock:: UTMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[i] += conj(L[n++]) * rhs[k];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: UMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[i] += U[n++] * rhs[k];
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: UTMult(Num* lhs, const Num* rhs, int nComp) const
{
    int i,k,n=1;
    FORNCOMP(i) FORNCOMP(k)  lhs[k] += conj(U[n++]) * rhs[i];
}
//-------------------------------------------------------------------------

Num& NeighbourBlock:: operator()(int comp1, int comp2, int nComp) 
{
    return L[nComp*(comp1-1) + comp2]; 
}
//-------------------------------------------------------------------------

Num& NeighbourBlock:: operator()(int comp1, int comp2, int nComp, 
				 Bool /*transpose*/) 
{
    return L[nComp*(comp2-1) + comp1]; 
}
//-------------------------------------------------------------------------

Num& AsymNeighbourBlock:: operator()(int comp1, int comp2, int nComp,
				     Bool /*transpose*/) 
{
    return U[nComp*(comp2-1) + comp1]; 
}
//-------------------------------------------------------------------------

void BlockNode:: setDataSpace(FixedSizeAllocator& dataAlloc, int nComp)
{
    if (D)    { cout << "\n*** NeighbourBlock: D != 0 \n"; cout.flush(); abort(); }
    if (DInv) { cout << "\n*** NeighbourBlock: DInv != 0 \n"; cout.flush(); abort(); }

    D    = (Num*) dataAlloc.Get();
    DInv = (Num*) dataAlloc.Get();
    --D;
    --DInv;

    int i;
    const int size = nComp*nComp;
    for (i=1; i<=size; ++i) D[i] = DInv[i] = 0.0; 
}
//-------------------------------------------------------------------------


void BlockNode:: removeDataSpace(FixedSizeAllocator& dataAlloc)
{
    ++D;    dataAlloc.Return(D);     D    = 0;
    ++DInv; dataAlloc.Return(DInv);  DInv = 0;
}
//-------------------------------------------------------------------------


void NeighbourBlock:: setDataSpace(FixedSizeAllocator& dataAlloc, int nComp)
{
    if (L) { cout << "\n*** NeighbourBlock: L != 0 \n"; cout.flush(); abort(); }
    L = (Num*) dataAlloc.Get();
    --L;

    int i;
    const int size = nComp*nComp;
    for (i=1; i<=size; ++i) L[i] = 0.0; 
}
//-------------------------------------------------------------------------


void NeighbourBlock:: removeDataSpace(FixedSizeAllocator& dataAlloc)
{
    ++L;  dataAlloc.Return(L);  L = 0;
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: setDataSpace(FixedSizeAllocator& dataAlloc, int nComp)
{
    if (L) { cout << "\n*** NeighbourBlock: L != 0 \n"; cout.flush(); abort(); }
    if (U) { cout << "\n*** NeighbourBlock: U != 0 \n"; cout.flush(); abort(); }

    L = (Num*) dataAlloc.Get();
    U = (Num*) dataAlloc.Get();
    --L;
    --U;

    int i;
    const int size = nComp*nComp;
    for (i=1; i<=size; ++i) L[i] = U[i] = 0.0; 
}
//-------------------------------------------------------------------------

void AsymNeighbourBlock:: removeDataSpace(FixedSizeAllocator& dataAlloc)
{
    ++L;  dataAlloc.Return(L);  L = 0;
    ++U;  dataAlloc.Return(U);  U = 0;
}