/home/nutarojj/Code/thyme-2.1.1/src/AdmittanceNetwork_SUPERLU.cpp

#include "AdmittanceNetwork_SUPERLU.h"
#include <iostream>
using namespace std;

AdmittanceNetwork_SUPERLU::AdmittanceNetwork_SUPERLU(int N, bool init_eye):
        N(N), // Remember the size of the matrix
        Nalloc(N)
{
        // Allocate the permutation matrices
        if ( !(perm_r = intMalloc(N)) ) ABORT("Malloc failed for perm_r[]");
        if ( !(perm_c = intMalloc(N)) ) ABORT("Malloc failed for perm_c[]");
        // Make sure all matrices are computed at
        // first call to solve_for_voltage
        dirty_val = dirty_struct = true;
        // Allocate space for and initialize the b vector
        if ( !(rhs = doublecomplexMalloc(N)) ) ABORT("Malloc failed for rhs[]");
        zCreate_Dense_Matrix(&B,N,1,rhs,N,SLU_DN,SLU_Z,SLU_GE);
        // Allocate space for and initialize the A matrix
        doublecomplex *a;
        int *arow, *acol;
        if ( !(a = doublecomplexMalloc(N)) ) ABORT("Malloc failed for a[]");
        if ( !(arow = intMalloc(N)) ) ABORT("Malloc failed for arow[]");
        if ( !(acol = intMalloc(N+1)) ) ABORT("Malloc failed for acol[]");
        for (int i = 0; i < N; i++)
        {
                // Initialize the diagonals with 1+j0 or 0+j0
                a[i].r = 1.0*init_eye; a[i].i = 0.0;
                arow[i] = acol[i] = i;
        }
        acol[N] = N;
        zCreate_CompCol_Matrix(&A,N,N,N,a,arow,acol,SLU_NC,SLU_Z,SLU_GE);
        // L and U are created when the matrix is factored
        L.Store = U.Store = NULL;
        // Options for the LU solver
        set_default_options(&options);
}

AdmittanceNetwork_SUPERLU::~AdmittanceNetwork_SUPERLU()
{
        SUPERLU_FREE(rhs);
        SUPERLU_FREE(perm_r);
        SUPERLU_FREE(perm_c);
        Destroy_CompCol_Matrix(&A);
        Destroy_SuperMatrix_Store(&B);
        if (L.Store != NULL) // Do these exist?
        {
                Destroy_SuperNode_Matrix(&L);
                Destroy_CompCol_Matrix(&U);
        }
}

void AdmittanceNetwork_SUPERLU::set(int i, int j, Complex z)
{
        // Get the data holding the matrix entries
        NCformat *Astore = (NCformat*)(A.Store);
        doublecomplex *data = (doublecomplex*)(Astore->nzval);
        // Look for an existing entry
        const int col_end = Astore->colptr[j+1];
        for (int col = Astore->colptr[j]; col < col_end; col++)
        {
                // If the entry is found
                if (Astore->rowind[col] == i)
                {
                        // Was the value changed?
                        bool new_dirty_val = (data[col].r != real(z) ||
                                        data[col].i != imag(z));
                        if (!new_dirty_val) return;
                        dirty_val = true;
                        // If so, assign a new value
                        data[col].r = real(z);
                        data[col].i = imag(z);
                        // If this value is zero, then erase the entry
                        if (data[col].i == 0.0 && data[col].r == 0.0)
                        {
                                // Shift value and row indices to the left
                                for (int p = col+1; p < Astore->nnz; p++)
                                {
                                        data[p-1].r = data[p].r;
                                        data[p-1].i = data[p].i;
                                        Astore->rowind[p-1] = Astore->rowind[p];
                                }
                                // Shift the columns left by one to account
                                // for the deleted entry
                                for (int p = j+1; p < N; p++)
                                {
                                        Astore->colptr[p]--;
                                }
                                // Change the entry count and column end
                                Astore->colptr[N] = --(Astore->nnz);
                                // Note the change in structure
                                dirty_struct = true;
                        }
                        return;
                }
        }
        // Arriving here means that no existing entry was found and one must
        // be created, or if the new value is zero then just return
        if (real(z) == 0.0 && imag(z) == 0.0) return;
        // Add the new entry
        Astore->colptr[N] = ++(Astore->nnz);
        // Note the change
        dirty_val = dirty_struct = true;
        // Allocate additional space
        if (Astore->nnz >= Nalloc)
        {
                Nalloc += 1;
                Astore->rowind =
                        (int*)realloc(Astore->rowind,Nalloc*sizeof(int));
                Astore->nzval = realloc(Astore->nzval,
                                        Nalloc*sizeof(doublecomplex));
                data = (doublecomplex*)(Astore->nzval);
        }
        // Shift all of the stored values to the right
        for (int p = Astore->nnz-1; p > col_end; p--)
        {
                Astore->rowind[p] = Astore->rowind[p-1];
                data[p] = data[p-1];
        }
        // Put the new value into its position
        data[col_end].r = real(z);
        data[col_end].i = imag(z);
        Astore->rowind[col_end] = i;
        // Advance all of the column starting points after the changed column by one
        for (int p = j+1; p < N; p++) Astore->colptr[p]++;
}

Complex AdmittanceNetwork_SUPERLU::get(int i, int j) const
{
        Complex result(0.0,0.0);
        NCformat* Astore = (NCformat*)(A.Store);
        doublecomplex *data = (doublecomplex*)(Astore->nzval);
        int col = Astore->colptr[j];
        int col_end = Astore->colptr[j+1];
        for (; col < col_end; col++)
        {
                if (Astore->rowind[col] == i)
                {
                        result = Complex(data[col].r,data[col].i);
                        break;
                }
        }
        return result;
}

void AdmittanceNetwork_SUPERLU::add_line(int i, int j, Complex y, double turns, double phase_shift, int side)
{
        Complex a = polar(turns,phase_shift);
        set(side,side,get(side,side)+a*a*y);
        int off_side = i;
        if (side == off_side) off_side = j;
        set(off_side,off_side,get(off_side,off_side)+y);
        Complex tmp = get(i,j)-a*y;
        set(i,j,tmp);
        set(j,i,tmp);
}

void AdmittanceNetwork_SUPERLU::remove_line(int i, int j, Complex y, double turns, double phase_shift, int side)
{
        Complex a = polar(turns,phase_shift);
        set(side,side,get(side,side)-a*a*y);
        int off_side = i;
        if (side == off_side) off_side = j;
        set(off_side,off_side,get(off_side,off_side)-y);
        Complex tmp = get(i,j)+a*y;
        set(i,j,tmp);
        set(j,i,tmp);
}

void AdmittanceNetwork_SUPERLU::add_line(int i, int j, Complex y)
{
        set(i,i,get(i,i)+y);
        set(j,j,get(j,j)+y);
        // Can avoid a second call for j,i because the matrix
        // is symettric
        Complex tmp = get(i,j)-y; 
        set(i,j,tmp);
        set(j,i,tmp);
}

void AdmittanceNetwork_SUPERLU::remove_line(int i, int j, Complex y)
{
        add_line(i,j,-y);
}

void AdmittanceNetwork_SUPERLU::add_self(int i, Complex y)
{
        set(i,i,get(i,i)+y);
}

void AdmittanceNetwork_SUPERLU::remove_self(int i, Complex y)
{
        add_self(i,-y);
}

void AdmittanceNetwork_SUPERLU::solve_for_voltage(const Complex* I, Complex* V)
{
        // Copy current to the X matrix
        for (int i = 0; i < N; i++)
        {
                rhs[i].r = real(I[i]);
                rhs[i].i = imag(I[i]);
        }
        // Solve the equations
        int info;
        SuperLUStat_t stat;
        StatInit(&stat);
        // Matrix needs to be completely refactored
        if (dirty_struct || dirty_val)
        {
                options.Fact = DOFACT;
                if (U.Store != NULL)
                {
                        Destroy_CompCol_Matrix(&U);
                        Destroy_SuperNode_Matrix(&L);
                }
                zgssv(&options,&A,perm_c,perm_r,&L,&U,&B,&stat,&info);
        }
        // Matrix is already factored; just solve for x
        else
        {
                options.Fact = FACTORED;
                zgstrs(NOTRANS,&L,&U,perm_c,perm_r,&B,&stat,&info);
        }
        StatFree(&stat);
        dirty_val = dirty_struct = false;
        // Copy result to V
        for (int i = 0; i < N; i++)
        {
                V[i] = Complex(rhs[i].r,rhs[i].i);
        }
}

void AdmittanceNetwork_SUPERLU::solve_for_current(const Complex* V, Complex* I)
{
        // Zero I
        for (int p = 0; p < N; p++)
        {
                I[p] = Complex(0.0,0.0);
        }
        // Get the data holding the matrix entries
        NCformat *Astore = (NCformat*)(A.Store);
        doublecomplex *data = (doublecomplex*)(Astore->nzval);
        // Perform the multiplication
        int col = 0; // Current column
        for (int p = 0; p < Astore->nnz; p++)
        {
                // Move to the next column?
                if (Astore->colptr[col+1] == p) col++;
                // Update I
                Complex a(data[p].r,data[p].i);
                I[Astore->rowind[p]] += a*V[col];
        }
}

---