Sparse.rar_matlab稀疏矩阵_sparse_稀疏优化_稀疏矩阵_节点优化

function [LCO,Diag,URO]=LUFact(A,Cindx,ERP,ERPU,CindxUO,Nbus) %************************************************************************************************* %* LDU Factorization %************************************************************************************************* %* INPUT VARIABLES %* A: The admittence matrix %* Cindx:Column indices %* ERP: End of row pointer %* Nbus: number of buses %* ERPU: the end of row pointer of upper triangular matirx %* CindxUO: Ordered column indices of rows in upper triangular matrix %************************************************************************************************* %* OUTPUT VECTOR %* LCO:Lower Triangular Matrix %* Diag:Diagonal Matrix %* URO:Upper Triangular Matrix %************************************************************************************************* %* INTERNAL VARIABLES %* i,j: Counter %* t: Use to find element in Link %* t1: Use to temporarily store Link element (to deal with disjoint) %* URO: Elements of upper triangular matrix (stored by row) %* LCO: Elements of lower triangular matrix (stored by columns) %* Diag: Inversed elements of diagonal matrix %* Diag: Initial collumn pointers of L %* ExAcum:Expanded Accumulator used in update process %* Link: Self-referential linked list used to associate %* the disjoint sets of rows needed in the reduction of each row %* Rindex:Row index of current row being processed %************************************************************************************************* %************************************************************************** %* Initialize URO, LCO, ICPL, Diag, Link, ExAcum %************************************************************************** URO=zeros(1,ERPU(Nbus-1)); LCO=zeros(1,ERPU(Nbus-1)); ICPL=zeros(1,Nbus-1); Diag=zeros(1,Nbus); Link=zeros(1,Nbus); ExAcum=zeros(1,Nbus); %************************************************************************** %* Initialize ICPL from ERPU %************************************************************************** ICPL(1)=1; for i=2:Nbus-1 ICPL(i)=ERPU(i-1)+1; end %************************************************************************** %* Reduction from row 1 to Nbus-1 %************************************************************************** for i=1:Nbus-1 Rindex=i; %========================================================================= % Row 1 %========================================================================= %%************************************************************************ %%* Step 1: %%* Step 1.a & 1.b: Zero ExAcum using Link and CindxUO %%* Step 1.c: Copy elements from row Rindex of A to ExAcum %%************************************************************************* if Rindex==1 % for row 1 the pointer from 1 to ERPU(1) for j=1:ERPU(Rindex) % Zero ExAcum using CindxUO ExAcum(CindxUO(j))=0; end t=Rindex; % Zero ExAcum using Link while Link(t)~=0 ExAcum(Link(t))=0; t=Link(t); end for j=1:ERP(Rindex) % Copy elements from A to ExAcum ExAcum(Cindx(j))=A(j); end %%************************************************************************ %%* Step 2 %%************************************************************************* %Step 2.a RX=0; %Initialize RX %Step 2.b Search Link list while Link(Rindex)~=0 RX=Link(Rindex); % Step 2.d Multiply row RX (stored in URO) by ExAcum(RX) % and subtract from appropriate locations in ExAcum if RX==1 for j=1:ERPU(RX) ExAcum(CindxUO(j))=ExAcum(CindxUO(j))-URO(j)*ExAcum(RX); end else for j=ERPU(RX-1)+1:ERPU(RX) ExAcum(CindxUO(j))=ExAcum(CindxUO(j))-URO(j)*ExAcum(RX); end end %*************************************************************** % Step 2.e & f Multiply Diag(RX) by ExAcum(RX) % and storeed into L matrix by ICPL(RX) LCO(ICPL(RX))=ExAcum(RX)*Diag(RX); ICPL(RX)=ICPL(RX)+1; %************************************************************** %* Disjoint the row RX from Link list associate with Row Rindex %************************************************************** t=Link(Rindex); t1=Link(t); Link(t)=0; Link(Rindex)=t1; %************************************************************** %Step 2.g Update the Link list add RX to the Link associate with % Row CindxUO(ICPL(RX)) if ICPL(RX)<=ERPU(RX) t=CindxUO(ICPL(RX)); while Link(t)>0 t=Link(t); end Link(t)=RX; end end %************************************************************************** %* Step 3: Zero all elements in Link Associate with row Rindex %************************************************************************** t=Rindex; while Link(t)>0 t1=Link(t); Link(t)=0; t=t1; end %************************************************************************** %* Step 4: Invert diagonal element and stored in Diag(Rindex) %************************************************************************** Diag(Rindex)=1/ExAcum(Rindex); %************************************************************************** %* Step 5: Update URO by multiply the ExAcum(CindxUO(k)) with Diag(Rindex) %************************************************************************** for j=1:ERPU(Rindex) URO(j)=ExAcum(CindxUO(j))*Diag(Rindex); end %************************************************************************** %* Step 6: Add Rindex into Link list %************************************************************************** t=CindxUO(ICPL(Rindex)); while Link(t)>0 t=Link(t); end Link(t)=Rindex; %************************************************************************* %************************************************************************* %************************************************************************* %************************************************************************* %========================================================================= % Row > 1 %========================================================================= else %%************************************************************************ %%* Step 1: %%* Step 1.a & 1.b: Zero ExAcum using Link and CindxUO %%* Step 1.c: Copy elements from row Rindex of A to ExAcum %%************************************************************************* % Zero ExAcum using CindxUO for j=ERPU(Rindex-1)+1:ERPU(Rindex) % for row Rindex the pointer from ERPU(Rindex-1)+1 to ERPU(Rindex) ExAcum(CindxUO(j))=0; end % Zero ExAcum using Link t=Rindex; while Link(t)~=0 ExAcum(Link(t))=0; t=Link(t); end % Copy elements from A to ExAcum for j=ERP(Rindex-1)+1:ERP(Rindex) ExAcum(Cindx(j))=A(j); end %%************************************************************************ %%* Step 2 %%************************************************