00001 #include "AdmittanceNetwork_SUPERLU.h"
00002 #include <iostream>
00003 using namespace std;
00004
00005 AdmittanceNetwork_SUPERLU::AdmittanceNetwork_SUPERLU(int N, bool init_eye):
00006 N(N),
00007 Nalloc(N)
00008 {
00009
00010 if ( !(perm_r = intMalloc(N)) ) ABORT("Malloc failed for perm_r[]");
00011 if ( !(perm_c = intMalloc(N)) ) ABORT("Malloc failed for perm_c[]");
00012
00013
00014 dirty_val = dirty_struct = true;
00015
00016 if ( !(rhs = doublecomplexMalloc(N)) ) ABORT("Malloc failed for rhs[]");
00017 zCreate_Dense_Matrix(&B,N,1,rhs,N,SLU_DN,SLU_Z,SLU_GE);
00018
00019 doublecomplex *a;
00020 int *arow, *acol;
00021 if ( !(a = doublecomplexMalloc(N)) ) ABORT("Malloc failed for a[]");
00022 if ( !(arow = intMalloc(N)) ) ABORT("Malloc failed for arow[]");
00023 if ( !(acol = intMalloc(N+1)) ) ABORT("Malloc failed for acol[]");
00024 for (int i = 0; i < N; i++)
00025 {
00026
00027 a[i].r = 1.0*init_eye; a[i].i = 0.0;
00028 arow[i] = acol[i] = i;
00029 }
00030 acol[N] = N;
00031 zCreate_CompCol_Matrix(&A,N,N,N,a,arow,acol,SLU_NC,SLU_Z,SLU_GE);
00032
00033 L.Store = U.Store = NULL;
00034
00035 set_default_options(&options);
00036 }
00037
00038 AdmittanceNetwork_SUPERLU::~AdmittanceNetwork_SUPERLU()
00039 {
00040 SUPERLU_FREE(rhs);
00041 SUPERLU_FREE(perm_r);
00042 SUPERLU_FREE(perm_c);
00043 Destroy_CompCol_Matrix(&A);
00044 Destroy_SuperMatrix_Store(&B);
00045 if (L.Store != NULL)
00046 {
00047 Destroy_SuperNode_Matrix(&L);
00048 Destroy_CompCol_Matrix(&U);
00049 }
00050 }
00051
00052 void AdmittanceNetwork_SUPERLU::set(int i, int j, Complex z)
00053 {
00054
00055 NCformat *Astore = (NCformat*)(A.Store);
00056 doublecomplex *data = (doublecomplex*)(Astore->nzval);
00057
00058 const int col_end = Astore->colptr[j+1];
00059 for (int col = Astore->colptr[j]; col < col_end; col++)
00060 {
00061
00062 if (Astore->rowind[col] == i)
00063 {
00064
00065 dirty_val = (data[col].r != real(z) ||
00066 data[col].i != imag(z));
00067 if (!dirty_val) return;
00068
00069 data[col].r = real(z);
00070 data[col].i = imag(z);
00071
00072 if (data[col].i == 0.0 && data[col].r == 0.0)
00073 {
00074
00075 for (int p = col+1; p < Astore->nnz; p++)
00076 {
00077 data[p-1].r = data[p].r;
00078 data[p-1].i = data[p].i;
00079 Astore->rowind[p-1] = Astore->rowind[p];
00080 }
00081
00082
00083 for (int p = j+1; p < N; p++)
00084 {
00085 Astore->colptr[p]--;
00086 }
00087
00088 Astore->colptr[N] = --(Astore->nnz);
00089
00090 dirty_struct = true;
00091 }
00092 return;
00093 }
00094 }
00095
00096
00097 if (real(z) == 0.0 && imag(z) == 0.0) return;
00098
00099 Astore->colptr[N] = ++(Astore->nnz);
00100
00101 dirty_val = dirty_struct = true;
00102
00103 if (Astore->nnz >= Nalloc)
00104 {
00105 Nalloc += 1;
00106 Astore->rowind =
00107 (int*)realloc(Astore->rowind,Nalloc*sizeof(int));
00108 Astore->nzval = realloc(Astore->nzval,
00109 Nalloc*sizeof(doublecomplex));
00110 data = (doublecomplex*)(Astore->nzval);
00111 }
00112
00113 for (int p = Astore->nnz-1; p > col_end; p--)
00114 {
00115 Astore->rowind[p] = Astore->rowind[p-1];
00116 data[p] = data[p-1];
00117 }
00118
00119 data[col_end].r = real(z);
00120 data[col_end].i = imag(z);
00121 Astore->rowind[col_end] = i;
00122
00123 for (int p = j+1; p < N; p++) Astore->colptr[p]++;
00124 }
00125
00126 Complex AdmittanceNetwork_SUPERLU::get(int i, int j) const
00127 {
00128 Complex result(0.0,0.0);
00129 NCformat* Astore = (NCformat*)(A.Store);
00130 doublecomplex *data = (doublecomplex*)(Astore->nzval);
00131 int col = Astore->colptr[j];
00132 int col_end = Astore->colptr[j+1];
00133 for (; col < col_end; col++)
00134 {
00135 if (Astore->rowind[col] == i)
00136 {
00137 result = Complex(data[col].r,data[col].i);
00138 break;
00139 }
00140 }
00141 return result;
00142 }
00143
00144 void AdmittanceNetwork_SUPERLU::add_line(int i, int j, Complex y, double turns, double phase_shift, int side)
00145 {
00146 Complex a = polar(turns,phase_shift);
00147 set(side,side,get(side,side)+a*a*y);
00148 int off_side = i;
00149 if (side == off_side) off_side = j;
00150 set(off_side,off_side,get(off_side,off_side)+y);
00151 Complex tmp = get(i,j)-a*y;
00152 set(i,j,tmp);
00153 set(j,i,tmp);
00154 }
00155
00156 void AdmittanceNetwork_SUPERLU::remove_line(int i, int j, Complex y, double turns, double phase_shift, int side)
00157 {
00158 Complex a = polar(turns,phase_shift);
00159 set(side,side,get(side,side)-a*a*y);
00160 int off_side = i;
00161 if (side == off_side) off_side = j;
00162 set(off_side,off_side,get(off_side,off_side)-y);
00163 Complex tmp = get(i,j)+a*y;
00164 set(i,j,tmp);
00165 set(j,i,tmp);
00166 }
00167
00168 void AdmittanceNetwork_SUPERLU::add_line(int i, int j, Complex y)
00169 {
00170 set(i,i,get(i,i)+y);
00171 set(j,j,get(j,j)+y);
00172
00173
00174 Complex tmp = get(i,j)-y;
00175 set(i,j,tmp);
00176 set(j,i,tmp);
00177 }
00178
00179 void AdmittanceNetwork_SUPERLU::remove_line(int i, int j, Complex y)
00180 {
00181 add_line(i,j,-y);
00182 }
00183
00184 void AdmittanceNetwork_SUPERLU::add_self(int i, Complex y)
00185 {
00186 set(i,i,get(i,i)+y);
00187 }
00188
00189 void AdmittanceNetwork_SUPERLU::remove_self(int i, Complex y)
00190 {
00191 add_self(i,-y);
00192 }
00193
00194 void AdmittanceNetwork_SUPERLU::solve_for_voltage(const Complex* I, Complex* V)
00195 {
00196
00197 for (int i = 0; i < N; i++)
00198 {
00199 rhs[i].r = real(I[i]);
00200 rhs[i].i = imag(I[i]);
00201 }
00202
00203 int info;
00204 SuperLUStat_t stat;
00205 StatInit(&stat);
00206
00207 if (dirty_struct || dirty_val)
00208 {
00209 options.Fact = DOFACT;
00210 if (U.Store != NULL)
00211 {
00212 Destroy_CompCol_Matrix(&U);
00213 Destroy_SuperNode_Matrix(&L);
00214 }
00215 zgssv(&options,&A,perm_c,perm_r,&L,&U,&B,&stat,&info);
00216 }
00217
00218 else
00219 {
00220 options.Fact = FACTORED;
00221 zgstrs(NOTRANS,&L,&U,perm_c,perm_r,&B,&stat,&info);
00222 }
00223 StatFree(&stat);
00224 dirty_val = dirty_struct = false;
00225
00226 for (int i = 0; i < N; i++)
00227 {
00228 V[i] = Complex(rhs[i].r,rhs[i].i);
00229 }
00230 }
00231
00232 void AdmittanceNetwork_SUPERLU::solve_for_current(const Complex* V, Complex* I)
00233 {
00234
00235 for (int p = 0; p < N; p++)
00236 {
00237 I[p] = Complex(0.0,0.0);
00238 }
00239
00240 NCformat *Astore = (NCformat*)(A.Store);
00241 doublecomplex *data = (doublecomplex*)(Astore->nzval);
00242
00243 int col = 0;
00244 for (int p = 0; p < Astore->nnz; p++)
00245 {
00246
00247 if (Astore->colptr[col+1] == p) col++;
00248
00249 Complex a(data[p].r,data[p].i);
00250 I[Astore->rowind[p]] += a*V[col];
00251 }
00252 }
00253