CHAPTER 7: MATRIX FACTORIZATION AND OPTIMAL ORDERING

7.1 Triangularization of a Matrix
7.1 Triangularization of Matrix
7.2 Ordinary LU Factorization
7.3 LU-Dolittle Method (LUDM) for Sparse Matrix
7.4 Code Fragments on LUDM
7.5 Example Calculation
7.6 Optimal Ordering (OO)
7.7 OO Scheme_123 Computer Program (SCHEME123.C)

Triangularization is the process of applying Gauss-Jordan Reduction to solve the simultaneous linear equation of the form Ax=b. LU Factorization where L="Lower" and U="Upper" is synonymous to triangularization. Factorization simply means that when we multiply matrices L and U, we get matrix A. Given the LU, we can perform forward course and backward course on b to solve for x. This is needed in Loadflow calculation. With LU we can also get the inverse of A which is used in short circuit calculation. Factorization is used to solve large sparse matrix A. Sparse means that more A(i,j) elements are zeroes than are non-zeroes. For power systems, it can be just 2 per cent of the square matrix A is non-zero. There are many references on LU Factorization that exploits matrix sparsity. Ref [21] contains clear formula derivation.

7.2 Ordinary LU Factorization

We can derive the formula in LU factorization using a 4x4 size matrix. The problem can be formulated as:
Given matrix A, solve the elements in matrices L and U such that L*U = A,
where matrices L and U are factors of A.
Once the elements in L and U are solved, unknown vector x in the equation A*x = b can be determined by forward course and backward course using L and U . There are variations to the structure of L and U and other References call these LDU. We choose our form as shown below.

--- Eq. (1)

Performing matrix multiplication, each A(i,j) elements can be equated as product of a row of L and a column of U. Remember that matrix A is know quantity and matrices L and U are unknown factors of A.

Row=1. Equating elements on row=1 as shown by the dotted arrow. These elements are A(1,1), A(1,2), A(1,3), A(1,4):
A(1,1)=1*U(1,1); A(1,2)=1*U(1,2); A(1,3)=1*U(1,3);
        A(1,4)=1*U(1,4)
so U(1,1)=A(1,1)
    U(1,2)=A(1,2)
    U(1,3)=A(1,3)
    U(1,4)=A(1,4)
Col=1. Equating elements on col=1, these are A(2,1), A(3,1), A(4,1):
A(2,1)=L(2,1)*U(1,1); A(3,1)=L(3,1)*U(1,1);
        A(4,1)=L(4,1)*U(1,1)

At this stage, U(1,1) had been solved already. So
L(2,1)=A(2,1)/U(1,1)
L(3,1)=A(3,1)/U(1,1)
L(4,1)=A(4,1)/U(1,1)

Row=2. Equating elements A(2,2), A(2,3), A(2,4):
A(2,2)=L(2,1)*U(1,2) + 1*U(2,2)
so U(2,2)=A(2,2) - L(2,1)*U(1,2) where L(2,1) and U(1,2) had been previously solved.
Also,
A(2,3)=L(2,1)*U(1,3) + 1*U(2,3),

so U(2,3)=A(2,3) - L(2,1)*U(1,3), where U(1,3) was previously solved.
A(2,4)=L(2,1) - U(1,4) + 1*U(2,4),
so U(2,4)=A(2,4) - L(2,1)*U(1,4), where U(1,4) was previously solved.

Col=2. Equating elements A(3,2), A(4,2):
A(3,2)=L(3,1)*U(1,2) + L(3,2)*U(2,2)
so L(3,2)=(A(3,2) - L(3,1)*U(1,2))/U(2,2) where necessary right hand terms were previously solved.
A(4,2)=L(4,1)*U(1,2) + L(4,2)*U(2,2)
L(4,2)=(A(4,2) - L(4,1)*U(1,2))/U(2,2)
Row=3. Equating elements A(3,3), A(3,4):
A(3,3)=L(3,1)*U(1,3) + L(3,2)*U(2,3) + 1*U(3,3)
The unknown here is U(3,3).

U(3,3)=A(3,3) - L(3,1)*U(1,3) - L(3,2)*U(2,3)
A(3,4)=L(3,1)*U(1,4) + L(3,2)*U(2,4) + U(3,4)
U(3,4)=A(3,4) - L(3,1)*U(1,4) - L(3,2)*U(2,4)

Col=3. Equating elements A(4,3):

A(4,3)=L(4,1)*U(1,3) + L(4,2)*U(2,3) + L(4,3)*U(3,3)
L(4,3)=(A(4,3) - L(4,1)*U(1,3) - L(4,2)*U(2,3))/U(3,3)
For ordinary LU, we can generalize that for matrix of size N,
--- Eq. (2)
Equating elements A(4,4):

A(4,4)=L(4,1)*U(1,4) + L(4,2)*U(2,4)
+L(4,3)*U(3,4) + U(4,4)
U(4,4)=A(4,4) - L(4,1)*U(1,4) - L(4,2)*U(2,4)
- L(4,3)*U(3,4)
And also ,

--- Eq. (3)
By looking at Eqs (2) and (3), we can conclude that there is no need to store L,U and A as three matrices. Use only one matrix A. After the factorization, matrix A is destroyed and replaced with L and U. The diagonal of matrix L is not stored.
To do: Write a computer program to perform Ordinary LU Factorization using only one matrix.

7.3 LU-Dolittle Method (LUDM) for Sparse Matrix.
This factorization is applicable for sparse matrices with half-storage. After solving for row 1 and col 1, modify those inside the enclosed dotted square (or rectangle),
Loop 1:
A'(2,2) = A(2,2) - L(2,1)*U(1,2),
A'(2,3) = A(2,3) - L(2,1)*U(1,3),
A'(2,4) = A(2,4) - L(2,1)*U(1,4),
A'(3,2) = A(3,2) - L(3,1)*U(1,2),
A'(3,3) = A(3,3) - L(3,1)*U(1,3),
A'(3,4) = A(3,4) - L(3,1)*U(1,4),
A'(4,2) = A(4,2) - L(4,1)*U(1,2),
A'(4,3) = A(4,3) - L(4,1)*U(1,3),
A'(4,4) = A(4,4) - L(4,1)*U(1,4).

Where U(2,2) = A'(2,2), U(2,3) = A'(2,3), and U(2,4) = A'(2,4)
Also L(3,2) = A'(3,2)/U(2,2) and L(4,2) = A'(4,2)/U(2,2)

Loop 2 :

A"(3,3) = A'(3,3) - L(3,2)*U(2,3),
A"(3,4) = A'(3,4) - L(3,2)*U(2,4),
A"(4,3) = A'(4,3) - L(4,2)*U(2,3),
A"(4,4) = A'(4,4) - L(4,2)*U(2,4).

Where U(3,3) = A"(3,3) , U(3,4) = A"(3,4)
and L(4,3) = A"(4,3)/U(3,3).
We can check from Ordinary LU that,
U(3,3) = A"(3,3)
= A'(3,3) - L(3,2)*U(2,3)
= A(3,3) - L(3,1)*U(1,3) - L3,2*U2,3
which satisfies the previously derived value.

For symmetric matrix there is no need of saving matrix L. In program implementation we make a k-loop from 1 to N-1. As shown on Fig_6 on the k_th loop, we need to update A(i,j)'s inside the big red dotted square (for colored print). To solve a particular A(i,j), we need L(i,k) and U(k,j) pointed to by the small dotted red line. However L(i,k)=U(k,i)/U(k,k), so to save memory there is no need to store matrix L because its value can be derived from U itself. For half-storage, storing only U, then,
--- Eq. (4)

--- Eq. (5)

The original matrix [A] is destroyed. Resulting [U] is stored on [A] itself, so only one matrix is used. This is similar to Shipley inversion. By looking at Eq. (5) and Fig_7, we are now on loop=k and the rectangles represents rhs nodes of k. There are 4 here and by Eq. (5) we must modify 4 diagonal elements and 6 off-diagonals, ij, y1,y2,y3,y4, and y5 as we modify only the upper triangle. If location represented by ij is not in lifo then it is a fill-in. The purpose of Node Ordering is to minimize (not really absolute minimum) the total fill-ins introduced during the factorization.

7.4 Code Fragments on LUDM

code=1: LUDM factorization Full Matrix, Upper Triangle only 1: for (k=1;k< nb;k++) 2: for(i=k+1;i<= nb;i++) 3: for (j=i;j<= nb;j++) <- note j=i 4: Y[i][j]-=Y[k][i]*Y[k][j]/Y[k][k]; code=2: Full matrix, Simulates sparse code 1: for (k=1;k< nb;k++) 2: for(i=k+1;i<= nb;i++) 3: {Y[i][i]-=Y[k][i]*Y[k][i]/Y[k][k]; 4: for (j=i+1;j<= nb;j++) <- note j=i+1 5: Y[i][j]-=Y[k][i]*Y[k][j]/Y[k][k]; 6: } code=3: LUDM Factorization of Sparse Ybus,Compact Numbers,Without Ordering 1: void Factor(void) 2: {int kk,i,j,ij,ki,kj,I,J,EndBus[20],Adr[20],counter; 3: /*All matrices had been reset (=zero) at Ybus()*/ 4: for (kk=1;kk< nb;kk++) 5: {counter=Get_RHSnodes(kk,EndBus,Adr); //number of RHS elements 6: if (counter) /*perform finite change on upper triangle only*/ 7: for (i=1;i<= counter;i++) 8: {I=EndBus[i];ki=Adr[i]; 9: YPd[I]-=YP[ki]*YP[ki]/YPd[kk]; // Diagonal part 10: if (i< counter) 11: for (j=i+1;j<= counter;j++) // Off-diagonal part 12: {J=EndBus[j];kj=Adr[j]; 13: ij=BrnSrNsrt(I,J); //branch search-insert J 14: YP[ij]-=YP[kj]*YP[ki]/YPd[kk]; 15: } /*for j*/ 16: } /*for i*/ 17: } /*for kk*/ 18: } code=4: LUDM Factorization of Sparse Ybus,Non-compact Numbers,With Ordering 1: void Factor(void) 2: {int k,kk,i,j,ij,ki,kj,I,J,EndBus[20],Adr[20],counter; 3: for (k=1;k< nb;k++) 4: {kk=NewOrder[k]; 5: counter=Get_RHSnodes(kk,EndBus,Adr); 6: if (counter) //perform finite change on upper triangle only 7: for (i=1;i<= counter;i++) 8: {I=EndBus[i];ki=Adr[i]; 9: YPd[BusPtr[I]]-=YP[ki]*YP[ki]/YPd[kk]; //Diagonal part 10: if (i< counter) 11: for (j=i+1;j<= counter;j++) //Off-diagonal part 12: {J=EndBus[j];kj=Adr[j]; 13: ij=BrnSrNsrt(I,J); //branch search-insert I->J 14: YP[ij]-=YP[kj]*YP[ki]/YPd[kk]; 15: } /*for j*/ 16: } /*for i*/ 17: } /*for k*/ 18: }

Explanation of codes.

Both codes=1 and 2 are implemented in full matrix notation, but computation is done only on the upper trangle. In code=1 the line Y[i][j]-=Y[k][i]*Y[k][j]/Y[k][k]; should have been Yij = Yij - Yik*Ykj/Ykk if the lower part was also solved. But since only the upper triangle was solved and we do not store the L matrix, Yki is used instead of Yik because k is less than i and also k is less than j. Note that the third loop starts at j=i, while code=2 has j starts at j=i+1. In code=1, when i=k+1 and j=i (j=k+1) represents that the Yij being updated is the diagonal or Yii. It is thus equal to Yii-=Yki*Yki/Ykk. For code=2 and 3,the update of the diagonal must be separate from the update of the off diagonal, because the diagonal is saved in a different memory location, i.e. code=2 is a full matrix that simulates code=3. We can do this by having j=i+1. Before going to the j-loop, the diagonal Yii must have been updated first.

In code=3, the End nodes of k are stored in EndBus[] and its address is stored in Adr[] starting from 1 to counter. Only RHS nodes are stored with or without ordering scheme. The j-loop starts from j=i+1 similar to code=2. Prior to j-loop, the diagonal Yii which is stored in a separate vector, Ypd[I], where I is the end node. The off diagonal Yij is computed as; YP[ij]-=YP[kj]*YP[ki]/YPd[kk];

where kj and ki are addresses stored in Adr[], while ij is searched in lifo. If absent, this becomes a fill.

When i=counter, there is no more off diagonals to be solved. Thats why a test "if (iJ has to first use more=start[BusPtr[I]].

The call to "Get_RHSnodes(kk,EndBus,Adr)" for code=3 and code=4 are the same. "kk" in code=4 is also a compact number.

7.5 Example Calculation.

Using the previous Ybus data, we now perform
the LU Factorization in half-storage.
[A] here is the Ybus matrix. We follow the program flow.

k=1, Loop=1, Elements to modify are:
             Y(2,2), Y(2,3), Y(2,4), Y(3,3),Y(3,4), Y(4,4)
  i=2
    j=2
     Y(2,2) = Y(2,2)   -  Y(1,2)*Y(1,2)/Y(1,1)
            = 11.6667 - (-6.6667)*(-6.6667)/36.6667 = 11.6667 - 1.21213
            = 10.4545
    j=3
     Y(2,3) = Y(2,3) - Y(1,3)*Y(1,2)/Y(1,1)
            = 0      since Y(2,3) =0, and  Y(1,3) = 0
    j=4           
     Y(2,4) = Y(2,4) - Y(1,4)*Y(1,2)/Y(1,1)
            =    0       - (-10.0)*(-6.6667)/36.6667
            = -1.81818  
  this is a fill-in, since Y(2,4) = 0 previously.
 i=3
    j=3
     Y(3,3) = Y(3,3) - Y(1,3)*Y(1,3)/Y(1,1) 
            =      15   -  0*0/36.6667
            = 15
    j=4
     Y(3,4) = Y(3,4) - Y(1,3)*Y(1,4)/Y(1,1) 
            =    -10  -   0*(-10)/36.6667
            = -10
 i=4
    j=4
     Y(4,4) = Y(4,4) - Y(1,4)*Y(1,4)/Y(1,1) 
            =     20     - (-10)*(-10)/36.6667   = 20 - 2.7272
            = 17.2727
k=2, Loop=2, Elements to modify are: Y(3,3),Y(3,4), Y(4,4)
 i=3
    j=3
     Y(3,3) = Y(3,3) - Y(2,3)*Y(2,3)/Y(2,2) 
            = 15   -   0      = 15
    j=4     
     Y(3,4) = Y(3,4) - Y(2,3)*Y(2,4)/Y(2,2) 
            = -10 - 0          = -10
 i=4
    j=4
     Y(4,4) = Y(4,4) - Y(2,4)*Y(2,4)/Y(2,2) 
            = 17.2727 - (-1.81818)*(-1.81818)/10.4545
            = 17.2727 - 0.316206
            = 16.9565
k=3, Loop=3, Elements to modify are: Y(4,4)
 i=4
   j=4
    Y(4,4) = Y(4,4) - Y(3,4)*Y(3,4)/Y(3,3)
           = 16.9565 - (-10)*(-10)/15   = 16.9565 - 6.6667
           = 10.2899

Factored Ybus=

	1	2	3	4
1	36.6667	-6.6667	0	-10.0
2		10.4545	0	-1.81818 *
3			15.0	-10.0
4				10.2899

* where Y_2,4 is a fill-in.

7.6 Optimal Ordering.

When the matrix to be factored is sparse, the order in which rows are processed affects the number of fill-ins. If the order (or bus numbering) is shown as in Fig_8, there will be 4 fill-ins. If however, the node with most connections is numbered last, corresponding factorization will insert no fills as shown in Fig_9. Great savings in arithmetic operations and computer storage can be achieved if we can minimize the fill-ins introduced, particularly in repeat solutions as in Fast Decoupled Loadflow. Repeat solution means that the factored matrix remains the same for several iterations. Tinney-Walker [19] discusses optimal ordering known as schemes 1,2 and 3. The program implementation may not be exactly as described in Ref. 19.

Scheme 1
Number the rows according to the number of nonzero off-diagonal terms before the factorization. When fill-ins are added to the sparse matrix after factorization, the off-diagonal terms would increase. In this scheme the rows with only one off-diagonal term are numbered first sequentially from index 1 to bus count, those with two terms, etc., and those with the most terms, last. This scheme does not take into account any of the subsequent effects of the process. The only information needed is a list of the number of nonzero terms in each row (or off-diagonals) of the original matrix(i.e. unfactored). The second version uses a little technique to speed up ordering.
Scheme 2
Number the rows so that at each of the process the next row to be operated upon is the one with the fewest nonzero terms. If more than one row meets this criterion, select any one. This scheme requires insertion of fills in the process as the bus is ordered. Subtract one from NumOffDiag[ ] to all its rhs nodes and add one for nodes k-m (the fill-in). Increase valence counter until all the nodes are ordered.

Scheme 3
Number the rows so that at each step of the process the next row to be operated upon is the one that will introduce the fewest fill-in. We can speed up computer time by simulating factorization and process busses that have rhs nodes=valence count and order the busses starting from minimum fill. As the bus is ordered insert its fill and subtract one from NumOffDiag[ ] to all its rhs nodes and add one for nodes k-m (the fill-in). Repeat the loop since other NumOffDiag[ ] had changed until no nodes are ordered. This time increase valence counter until all the nodes are processed. Scheme-3 is implemented as,
1. Nodes with valence=1 are renumbered rightaway.
2. Compute table of fills for unrenumbered nodes, only for valence<=k.
3. The next node to be renumbered is the node with minimum fill, (say fill=1).
Renumber all nodes as determined in (2) starting from minimum fill, insert the fill.
4. Loop back to 2. This looping will take several times as other nodes have NumOffDiag[ ]
becoming =valens.
5. increase valence count and Loop back to 2 until all nodes are renumbered

In application programs such as loadflow and short circuit, the user must have an option to choose between schemes. However,in short circuit there is no such thing as repeat solution compared to loadflow.
Question: How do you know that a branch is a fill-in or an impedance?
Answer: A branch is a fill-in if its address is above 1->LINES. All fill-ins are appended after address LINES.
Observation: Scheme 2 and 3 simulates LU factorization. Fill-ins are inserted. So during actual LU factorization, there is no need to insert fil-ins,otherwise there is a bug in the program.
Further Research: Dynamic programming can be applied to scheme 3. The choice which node will be renumbered next should contribute minimum fill. All the paths taken should arrive at a global minimum.

7.7 OO Scheme_123 Computer Program (SCHEM123.C)

1: /* Node Ordering Using Tinney-Walker Schemes  file=schem123.c
   2:    Copyright 1996. All Rights Reserved. July 25,1996
   3:    Joesel Pante
   4:    Assumptions
   5:    1. contiguous bus numbering
   6:    2. no parallel lines
   7: scheme_0:  no scheme   
   8: scheme_1a: ordered according to number of off-diagonals only 
   9: scheme_1b: uses a technique to speed up ordering     
  10: scheme_2: according to rhs nodes
  11: scheme_3a: only one bus with minimum fill
  12: scheme_3b; several busses with minimum fill
  13: scheme_3c: all candidate busses at valence=k,according to fill_count
  14: */
  15: #define LMAX 200  //line data array size 
  16: #define BMAX 50  //bus data array size 
  17: #include <stdio.h>
  18: typedef enum {false,true,} Boolean;
  19: int OldBus[BMAX],NewBus[BMAX],NumOffDiag[BMAX],start[BMAX];
  20: int next[2*LMAX],bus[2*LMAX];  
  21: int mcount,nl,nb,bus_cntr,valens;
  22: /* NumOffDiag[i]= stores valence count, number of branches connected to
  23:    node=i       
  24:    start[ ]=pointer to first end node
  25:    next[ ]=pointer to next end node
  26:    bus[ ]=bus number of end node    
  27:    nl=line count, nb=bus count;
  28:    bus_cntr=bus counter for the ordered new bus
  29:    valens=valence counter
  30:    
  31: */
  32: FILE *infile;
  33: FILE *outfile;
  34:    
  35: void Simu_LU(void); 
  36: 
  37: void enter(void)
  38: {char *ch; //fprintf(outfile,"Enter>\n");
  39:  //gets(ch);
  40: }  
  41:    
  42: void PrintOrder(void)
  43: {int i;
  44:  fprintf(outfile,"i: NOffD OldBus NewBus\n");
  45:  for (i=1;i<=nb;i++)
  46:    fprintf(outfile,"%3d%6d%6d%6d\n",
  47:       i,NumOffDiag[i],OldBus[i],NewBus[i]);
  48: }            
  49: 
  50: void PrintLifo(void)
  51: {int i;
  52:  fprintf(outfile,"Number of Lines=%d Largest bus=%d\n",nl,nb);
  53:  fprintf(outfile,"i:  start  next  bus\n");
  54:  for (i=1;i<=nb;i++)
  55:    fprintf(outfile,"%3d%6d%6d%6d\n",i,start[i],next[i],bus[i]);
  56:  for (i=nb+1;i<=2*nl;i++)
  57:    fprintf(outfile,"%3d%12d%6d\n",i,next[i],bus[i]);      
  58: }       
  59: 
  60: void input (void)  //called by each scheme
  61: {int k,m;
  62:  mcount=0;//initialize branch counter here for each scheme/read 
  63:  nb=0;nl=0;bus_cntr=0;       //initialize   
  64:  for (k=1;k<BMAX;k++)          //initialize
  65:   {NumOffDiag[k]=0;start[k]=0;OldBus[k]=0;NewBus[k
  66:   ]=0;}    
  67:  infile=fopen("order.dat","r");
  68:  fprintf(outfile,"Count k  m\n");  
  69:  while (1)
  70:   {nl++;fscanf(infile,"%d %d",&k,&m);
  71:    fprintf(outfile,"%4d%4d%4d\n",nl,k,m);
  72:    if (k==999) break;
  73:    next[++mcount]=start[k];start[k]=mcount;bus[mcount]=m;        
  74:    next[++mcount]=start[m];start[m]=mcount;bus[mcount]=k;        
  75: //increase valence count
  76:    NumOffDiag[k]++; NumOffDiag[m]++;
  77: //get the largest bus
  78:    if (k>nb) nb=k;
  79:    if (m>nb) nb=m;
  80:   } //end while(1)       
  81:  nl=mcount/2;   //actual line count  
  82:  k=fclose(infile); 
  83:  enter();
  84: }
  85:  
  86: void scheme_0(void)         //without a scheme
  87: {int i; 
  88:  fprintf(outfile,"Scheme_0\n");
  89:  input(); 
  90:  PrintLifo(); 
  91:  for (i=1;i<=nb;i++)
  92:   {OldBus[i]=i; NewBus[i]=i;
  93:   }   
  94:  PrintOrder();
  95:  enter(); 
  96:  Simu_LU();
  97:  i=mcount/2-nl;
  98:  fprintf(outfile,"fill-in=%3d\n",i); 
  99:  
 100: }         
 101:   
 102: void scheme_1a(void)
 103: {int i;
 104:  fprintf(outfile,"Scheme_1a\n"); 
 105:  valens=1;//start with valence=1 
 106:  input(); 
 107:  PrintLifo(); 
 108:  while (bus_cntr<nb)
 109:   {fprintf(outfile,"valence count=%3d   bus_cntr=%3d\n",valens,bus_cntr);
 110:    enter();
 111:    for (i=1;i<=nb;i++)
 112: // only OldBus[ ]=0 are candidate for ordering
 113: // check if NumOffDiag[i] is the same as the current valence 
 114:      {if (NumOffDiag[i]==valens)         
 115:         {OldBus[i]=++bus_cntr;NewBus[bus_cntr]=i;
 116:         }
 117:      }
 118:    valens++;  //see next valence     
 119:    PrintOrder();
 120:    }  
 121:  Simu_LU();
 122:  i=mcount/2-nl;
 123:  fprintf(outfile,"fill-in=%3d\n",i); 
 124: }        
 125:   
 126: void scheme_1b(void)   
 127: {int i,k,noffd,MaxOffDiag=0; 
 128:  int bound[15];    //valence group counter
 129: //Determine maximum off diagonal count,MaxOffDiag,
 130: //we declared matrix size=15 for bound[], so MaxOffDiag<15   
 131:  fprintf(outfile,"Scheme_1b\n");
 132:  input(); 
 133:  PrintLifo(); 
 134:  for (i=1;i<=nb;i++)
 135:    if (NumOffDiag[i]>MaxOffDiag) MaxOffDiag=NumOffDiag[i];
 136:  fprintf(outfile,"Maximum Off Diag=%3d\n",MaxOffDiag); 
 137: //zero out valence group counter
 138:  for (i=1;i<=MaxOffDiag;i++)
 139:     bound[i]=0; 
 140: //Count busses in a valence group
 141:  for (i=1;i<=nb;i++)
 142:    bound[NumOffDiag[i]]++;
 143: //Print bound count 
 144:  fprintf(outfile,"Initial bound\n");
 145:  for (i=1;i<=MaxOffDiag;i++)
 146:     fprintf(outfile,"bound(%2d)=%2d\n",i,bound[i]);
 147: //Accumulate count
 148:  for (i=2;i<=MaxOffDiag;i++)
 149:     bound[i]+=bound[i-1]; 
 150: //Print bound count again
 151:  fprintf(outfile,"After accumulate bound-1\n");
 152:  for (i=1;i<=MaxOffDiag;i++)
 153:     fprintf(outfile,"bound(%2d)=%2d\n",i,bound[i]);
 154: //Perform Ordering Here, no need for variable nbus  
 155:  for (i=nb;i>=1;i--)
 156:    {noffd=NumOffDiag[i];     //noffd=number of off-diagonals
 157:     k=bound[noffd];
 158:     OldBus[i]=k;      //i=old bus number,k=new bus number, 
 159:     NewBus[k]=i;
 160:     bound[noffd]--; //decrease bus count per valence group
 161:    }  
 162: //Print bound count after the order 
 163:  fprintf(outfile,"After ordering,bound\n");
 164:  for (i=1;i<=MaxOffDiag;i++)
 165:     fprintf(outfile,"bound(%2d)=%2d\n",i,bound[i]);  
 166: //
 167:  Simu_LU();
 168:  i=mcount/2-nl;
 169:  fprintf(outfile,"fill-in=%3d\n",i);     
 170: }    
 171: 
 172: void Simu_LU(void) //called by scheme 1/2
 173: {int ip,rhs_cntr,i,j,k,m,L,rhs_bus[15]; 
 174:  Boolean found;
 175:  for (ip=1;ip<=nb;ip++)
 176:   {m=NewBus[ip];     //
 177:    L=start[m];rhs_cntr=0; 
 178:    fprintf(outfile,"ip=%3d Bus=%3d\n",ip,m);
 179:    while (L)
 180:     {k=bus[L];
 181:      if (OldBus[m]<OldBus[k])   //m is lhs of k
 182:        rhs_bus[++rhs_cntr]=k;
 183:      L=next[L];
 184:     }
 185: //no insertion if rhs=1 only
 186:    if (rhs_cntr>1)
 187:     {for (i=1;i<rhs_cntr;i++)
 188:        {k=rhs_bus[i]; 
 189:            for (j=i+1;j<=rhs_cntr;j++)
 190:             {m=rhs_bus[j];
 191:              fprintf(outfile,"   searching k=%3d ,m=%3d\n",k,m);
 192: //search if bus k-m is in lifo
 193:              L=start[k];found=false;
 194:              while (L && !found)
 195:               {found=(m==bus[L]);L=next[L];}
 196: //insert if end of thread and not in lifo  
 197: //you can not insert k-m if k and/or m had been already ordered            
 198:              if (!found)                    
 199:               {next[++mcount]=start[k];start[k]=mcount;bus[mcount]=m;
 200:                next[++mcount]=start[m];start[m]=mcount;bus[mcount]=k;
 201:                fprintf(outfile,"**mcount=%3d inserts k=%3d m=%3d\n",mcount,k,m); 
 202:               // enter();
 203:               }
 204:             }//next j
 205:          }//next i  
 206:         }//if rhs     
 207:   }//next ip
 208: }    
 209:             
 210: void scheme1(void)//called by scheme_2/3 only
 211: {Boolean process=true;
 212:  int i,L,k;
 213: //processing (valens)=1  
 214:  fprintf(outfile,"ordering for valence=1\n");
 215:  valens=1;
 216:  while (process && bus_cntr<nb)
 217:   {process=false;
 218:    for (i=1;i<=nb;i++)     
 219:     {if (NumOffDiag[i]<=valens && !OldBus[i])  //check if candidate
 220:       {OldBus[i]=++bus_cntr;NewBus[bus_cntr]=i; //renumber
 221:        fprintf(outfile,"Ordered Oldbus=%2d Newbus=%2d\n",i,bus_cntr);
 222:        process=true; 
 223: //decrease valence count to all nodes adjacent to bus=i 
 224:        L=start[i];
 225:        while (L)
 226:         {k=bus[L];
 227:          if (!OldBus[k])  //check if i is lhs of k
 228:            {NumOffDiag[k]--;
 229:             fprintf(outfile,"-1 for bus=%2d NOffD=%2d\n",k,NumOffDiag[k]);
 230:            }                                       
 231:          L=next[L];
 232:         }
 233:        }//if renumber took place
 234:      }//for i       
 235:   } //while       
 236: }         
 237:  
 238: void insert(int Bus) 
 239: {int L,k,m,i,j,rhs_cntr,rhs_bus[20];
 240:  Boolean found; 
 241:  //initialize rhs counter
 242:  rhs_cntr=0;
 243:  L=start[Bus];
 244:  while (L)
 245:   {k=bus[L];
 246:    if (!OldBus[k])  //only unordered rhs nodes
 247:      {NumOffDiag[k]--;rhs_bus[++rhs_cntr]=k;
 248:       fprintf(outfile,"-1 for bus=%2d NOffD=%2d\n",k,NumOffDiag[k]);
 249:      }  //decrease NumOffDiag[k] because ip is lhs of k
 250:    L=next[L];
 251:   }       
 252: //No insertion if rhs=1 only
 253:  if (rhs_cntr>1)
 254:    {for (i=1;i<rhs_cntr;i++)
 255:      {k=rhs_bus[i]; 
 256:       for (j=i+1;j<=rhs_cntr;j++)
 257:        {m=rhs_bus[j];
 258: //search if bus k-m is in lifo
 259:         L=start[k]; found=false;
 260:         while (L && !found)
 261:           {found=(m==bus[L]);L=next[L];}
 262: //insert if end of thread and not in lifo  
 263: //you can not insert k-m if k and/or m had been already ordered            
 264:            if (!found)                    
 265:              {next[++mcount]=start[k];start[k]=mcount;bus[mcount]=m;
 266:               next[++mcount]=start[m];start[m]=mcount;bus[mcount]=k;
 267:               NumOffDiag[k]++;NumOffDiag[m]++; 
 268:               fprintf(outfile,"**mcount=%3d inserts k=%3d m=%3d\n",mcount,k,m);  
 269:               fprintf(outfile,"+1 for bus=%2d NOffD=%2d\n",k,NumOffDiag[k]); 
 270:               fprintf(outfile,"+1 for bus=%2d NOffD=%2d\n",m,NumOffDiag[m]);                
 271:              }
 272:            }//next j
 273:         }//next i  
 274:     }//if rhs 
 275: }
 276:         
 277: void scheme_2(void)
 278: {int i,ip; 
 279:  fprintf(outfile,"Scheme_2\n");
 280:  input(); 
 281:  scheme1();   //for valens=1  
 282:  enter();
 283:  PrintOrder();   
 284:  enter();
 285:  //simulate LU factorization
 286:  //initialize valens=2  
 287:  valens=2;
 288:  while (bus_cntr<nb)
 289:   {fprintf(outfile,"processing valence=%2d\n",valens);
 290:    for (ip=1;ip<=nb;ip++)
 291:     {fprintf(outfile,"ip=>%2d\n",ip);  
 292:      if (!OldBus[ip] && NumOffDiag[ip]<=valens) //order it
 293:       {OldBus[ip]=++bus_cntr;NewBus[bus_cntr]=ip;
 294:        fprintf(outfile,"Ordered Oldbus=%2d Newbus=%2d\n",ip,bus_cntr);
 295: //decrease valence count for rhs nodes
 296: //insert fill-in and increase valence
 297:      insert(ip);
 298:     }//if OldBus
 299:   } //next ip  
 300:   //enter();
 301:   PrintOrder();
 302: //increase valence for next pass                   
 303:   valens++;
 304:  } //while  
 305:  i=mcount/2-nl;
 306:  fprintf(outfile,"fill-in=%3d\n",i);  
 307:  i=fclose(outfile);
 308: }
 309: 
 310: void GetCandBus(int *CndtKtr, int *fill_min,int *fill_max,
 311:                 int fln_ctr[1],int cndt_bus[1])
 312: {int ip,j,i,L,m,k,rhs_cntr,rhs_bus[20];
 313:  Boolean found;
 314:  *CndtKtr=0;  //initialize candidate bus counter
 315:    for (ip=1;ip<=nb;ip++)
 316:     {if (!OldBus[ip] && NumOffDiag[ip]<=valens)
 317: //Candidate busses are those still unordered and has <=valence count     
 318:       {cndt_bus[++*CndtKtr]=ip;fln_ctr[*CndtKtr]=0; //initialize fill-count
 319: //Get number of rhs busses and simulate LUfac to determine fills
 320:        rhs_cntr=0;
 321:        L=start[ip];
 322:        while (L)
 323:         {k=bus[L];
 324:          if (!OldBus[k])  //only unordered rhs nodes are included
 325:            rhs_bus[++rhs_cntr]=k;
 326:            L=next[L];
 327:         }       
 328: //No insertion if rhs=1 only
 329:        if (rhs_cntr>1)
 330:         {for (i=1;i<rhs_cntr;i++)
 331:           {k=rhs_bus[i]; 
 332: //         fprintf(outfile,"i=%3d ,k=%3d\n",i,k);
 333:            for (j=i+1;j<=rhs_cntr;j++)
 334:             {m=rhs_bus[j];
 335: //             fprintf(outfile,"j=%3d ,m=%3d\n",j,m);
 336: //search if bus k-m is in lifo
 337:              L=start[k]; found=false;
 338:              while (L && !found)
 339:               {found=(m==bus[L]);L=next[L];}
 340: //fill-in if end of thread and not in lifo  
 341: //you can not insert k-m if k and/or m had been already ordered            
 342:              if (!found)                    
 343:                ++fln_ctr[*CndtKtr]; 
 344:             }//next j
 345:          }//next i  
 346:         }//if rhs 
 347:     }//if OldBus
 348:   } //next ip 
 349: //Print Candidate busses and fills
 350:   if (*CndtKtr)
 351:    {fprintf(outfile,"Candidates for renumbering are:\n");
 352:     for (i=1;i<=*CndtKtr;i++)
 353:      fprintf(outfile,"i=%3d bus=%3d fills=%3d\n",i,cndt_bus[i],fln_ctr[i]);
 354:     enter();
 355: //search for minimum/max fills
 356:     *fill_min=99;*fill_max=0;
 357:     for (i=1;i<=*CndtKtr;i++)
 358:      {if (fln_ctr[i]<*fill_min) *fill_min=fln_ctr[i];   
 359:       if (fln_ctr[i]>*fill_max) *fill_max=fln_ctr[i];
 360:      } 
 361:    }//if (*CndtKtr)     
 362: }
 363:         
 364: void scheme_3(int skim) //skim=1,a/=2,b/=3,c
 365: {int i,j,k,fill_min,fill_max;
 366:  int fln_ctr[50],CndtKtr,cndt_bus[50]; 
 367: //fln_ctr[]=fill-in counter
 368: //cndt_bus[]=candidate bus ,CndtKtr=count   
 369:  fprintf(outfile,"Scheme_3\n");
 370:  input(); 
 371:  scheme1();   //for valens=1
 372:  PrintOrder();   
 373:  enter();
 374:  //simulate LU factorization
 375:  //initialize valens=2  
 376:  valens=2;
 377:  while (bus_cntr<nb)
 378:   {fprintf(outfile,"processing valence=%2d\n",valens); 
 379: //CndtKtr is a switch, increase valens count if no more bus to process
 380:    do
 381:     {GetCandBus(&CndtKtr,&fill_min,&fill_max,fln_ctr,cndt_bus);
 382: //renumber according to 3a 3b 3c    
 383:      if (CndtKtr)  
 384:       {switch (skim)  
 385:        { case 1:       //order only one bus, the first with min fil
 386:          {for (i=1;i<=CndtKtr;i++)
 387:           if (fln_ctr[i]==fill_min) break;
 388:           j=cndt_bus[i];  
 389:           OldBus[j]=++bus_cntr;NewBus[bus_cntr]=j; //renumber
 390:           fprintf(outfile,"next ordered bus OldBus[%3d]=%3d\n",j,bus_cntr);
 391: //perform actual insertion 
 392:           insert(j);
 393:           enter();        
 394:           break;
 395:          }
 396:         case 2:    //order several busses with minimum fill
 397:         {for (i=1;i<=CndtKtr;i++)
 398:           if (fln_ctr[i]==fill_min)
 399:             {j=cndt_bus[i];  
 400: //Test for valens count because it will change as others are 
 401: //ordered ahead, inserting fills. Others are subtracted too if lhs
 402:              if (NumOffDiag[j]<=valens) 
 403:               {OldBus[j]=++bus_cntr;NewBus[bus_cntr]=j; //renumber
 404:                fprintf(outfile,"next ordered bus OldBus[%3d]=%3d\n",j,bus_cntr);
 405: //perform actual insertion 
 406:                insert(j);
 407:                enter();        
 408:               }//if NumOff
 409:             } 
 410:          break; 
 411:          }//case 1
 412:        
 413:         case 3:       //order all candidates
 414:          {for (k=fill_min;k<=fill_max;k++)
 415:             for (i=1;i<=CndtKtr;i++)
 416:             if (fln_ctr[i]==k)
 417:              {j=cndt_bus[i];  
 418: //Test for valens count because it will change as others are 
 419: //ordered ahead, inserting fills. Others are subtracted too if lhs
 420:               if (NumOffDiag[j]<=valens) 
 421:                {OldBus[j]=++bus_cntr;NewBus[bus_cntr]=j; //renumber
 422:                 fprintf(outfile,"next ordered bus OldBus[%3d]=%3d\n",j,bus_cntr);
 423: //perform actual insertion 
 424:                 insert(j);
 425:                // enter();        
 426:                }//if NumOff
 427:              } //if
 428:             }//for i
 429:            break;
 430:           }//case 3      
 431:         }//switch
 432:   } while (CndtKtr);
 433: //---next loop for valens=k but different fills 
 434: //or others whose valence decreased as nodes are ordered (lhs nodes)
 435: //increase valence for next pass  
 436:   valens++;fprintf(outfile,"increase valens counter=%2d\n",valens);
 437:   enter();
 438:  } //while (bus_cntr<nb) 
 439:  PrintOrder();
 440:  i=mcount/2-nl;
 441:  fprintf(outfile,"fill-in=%3d\n",i);  
 442:  i=fclose(outfile);
 443: }
 444: 
 445: void main(void)
 446: {int i;
 447:  outfile=fopen("output.txt","w"); 
 448:  //scheme_0();
 449:  enter();   
 450: // scheme_1a(); 
 451:  enter();      
 452: // scheme_1b(); 
 453:  enter();   
 454:  //scheme_2();
 455:  enter();
 456:  //scheme_3(1);
 457: // enter();
 458: // scheme_3(2);
 459: // enter();
 460:  scheme_3(3);
 461:  enter(); 
 462:  i=fclose(outfile);
 463:  }

HOME
CHAPTER ONE: INTRODUCTION
CHAPTER TWO: PROGRAMMING TOPICS
CHAPTER THREE: MATRIX INVERSION ROUTINE
CHAPTER FOUR: STEP-BY-STEP ALGORITHMS
CHAPTER FOUR: Part 2
CHAPTER FIVE: DIRECT METHODS IN Ybus FORMATION
CHAPTER SIX: SPARSITY
CHAPTER SEVEN: MATRIX FACTORIZATION AND OPTIMAL ORDERING
CHAPTER EIGHT: SPARSE Zbus FROM FACTORED SPARSE Ybus
CHAPTER NINE: FAULT CALCULATIONS
REFERENCES
Copyright � 1997 by Sparsity Technology. All rights reserved.

i n n o v a t e o r s t a g n a t e