00001 #include "AdmittanceNetwork_LAPACK.h"
00002 #include <iostream>
00003 #include <cassert>
00004 using namespace std;
00005
00006
00007 extern "C"
00008 {
00009 void zgetrs_(char*,int*,int*,double*,int*,int*,double*,int*,int*);
00010 void zgetrf_(int*,int*,double*,int*,int*,int*);
00011 };
00012
00013 AdmittanceNetwork_LAPACK::AdmittanceNetwork_LAPACK(int N):
00014 N(N),
00015 Y(new double[2*N*N]),
00016 LUY(new double[2*N*N]),
00017 dirty(false),
00018 iv_lapack(new double[2*N]),
00019 ipiv_lapack(new int[N])
00020 {
00021 for (int i = 0; i < 2*N*N; i++)
00022 Y[i] = 0.0;
00023 }
00024
00025 void AdmittanceNetwork_LAPACK::set(int i, int j, Complex y)
00026 {
00027 dirty = true;
00028 int real_part = i*2*N+2*j;
00029 int imag_part = real_part+1;
00030 Y[real_part] = real(y);
00031 Y[imag_part] = imag(y);
00032 }
00033
00034 Complex AdmittanceNetwork_LAPACK::get(int i, int j) const
00035 {
00036 int real_part = i*2*N+2*j;
00037 int imag_part = real_part+1;
00038 return Complex(Y[real_part],Y[imag_part]);
00039 }
00040
00041 void AdmittanceNetwork_LAPACK::add_line(int i, int j, Complex y)
00042 {
00043 set(i,i,get(i,i)+y);
00044 set(j,j,get(j,j)+y);
00045 set(i,j,get(i,j)-y);
00046 set(j,i,get(j,i)-y);
00047 }
00048
00049 void AdmittanceNetwork_LAPACK::add_line(int i, int j, Complex y, double turns, double phase_shift, int side)
00050 {
00051 Complex a = polar(turns,phase_shift);
00052 set(side,side,get(side,side)+a*a*y);
00053 int off_side = i;
00054 if (side == off_side) off_side = j;
00055 set(off_side,off_side,get(off_side,off_side)+y);
00056 Complex tmp = get(i,j)-a*y;
00057 set(i,j,tmp);
00058 set(j,i,tmp);
00059 }
00060
00061 void AdmittanceNetwork_LAPACK::remove_line(int i, int j, Complex y, double turns, double phase_shift, int side)
00062 {
00063 Complex a = polar(turns,phase_shift);
00064 set(side,side,get(side,side)-a*a*y);
00065 int off_side = i;
00066 if (side == off_side) off_side = j;
00067 set(off_side,off_side,get(off_side,off_side)-y);
00068 Complex tmp = get(i,j)+a*y;
00069 set(i,j,tmp);
00070 set(j,i,tmp);
00071 }
00072
00073 void AdmittanceNetwork_LAPACK::add_self(int i, Complex y)
00074 {
00075 set(i,i,get(i,i)+y);
00076 }
00077
00078 void AdmittanceNetwork_LAPACK::solve_for_current(const Complex* V, Complex* I)
00079 {
00080
00081 for (int i = 0; i < N; i++)
00082 {
00083 I[i] = Complex(0.0,0.0);
00084 for (int j = 0; j < N; j++)
00085 {
00086 I[i] += get(j,i)*V[j];
00087 }
00088 }
00089 }
00090
00091 void AdmittanceNetwork_LAPACK::solve_for_voltage(const Complex* I, Complex* V)
00092 {
00093 int i,info,SIZE=N,lda=N,ldb=N,nrhs=1;
00094 char type = 'N';
00095
00096 if (dirty)
00097 {
00098 dirty = false;
00099 for (int i = 0; i < 2*N*N; i++) LUY[i] = Y[i];
00100 zgetrf_(&SIZE,&SIZE,LUY,&lda,ipiv_lapack,&info);
00101 if (info != 0)
00102 {
00103 cerr << "zgetrf_ returned " << info << endl;
00104 cerr.flush();
00105 }
00106 assert(info == 0);
00107 }
00108
00109 for (i = 0; i < N; i++)
00110 {
00111 iv_lapack[2*i] = real(I[i]);
00112 iv_lapack[2*i+1] = imag(I[i]);
00113 }
00114
00115 zgetrs_(&type,&SIZE,&nrhs,LUY,&lda,ipiv_lapack,iv_lapack,&ldb,&info);
00116 if (info != 0)
00117 {
00118 cerr << "zgetrs_ returned " << info << endl;
00119 cerr.flush();
00120 }
00121 assert(info == 0);
00122
00123 for (i = 0; i < N; i++)
00124 {
00125 V[i] = Complex(iv_lapack[2*i],iv_lapack[2*i+1]);
00126 }
00127 }
00128
00129 AdmittanceNetwork_LAPACK::~AdmittanceNetwork_LAPACK()
00130 {
00131 delete [] Y;
00132 delete [] LUY;
00133 delete [] iv_lapack;
00134 delete [] ipiv_lapack;
00135 }
00136