Thread: LUP Decomposition (simultaneous equations)

  1. #1
    Registered User Markallen85's Avatar
    Join Date
    Nov 2002
    Posts
    53

    LUP Decomposition (simultaneous equations)

    I'm in a situation where I have to be able to solve an unknown number of simultaneous equations, there could be from two variables to a near unlimited number (practical limits will probably be in the 10-15 area, but this program contains a lot of random events, so ther could potentially be much more).

    The equations are all simple linear equations, and there will always be N equations, where N is the number of variables to solve. So the situation is not overly complex, however I need to do this fast (possibly millions of times in s single execution) so I want to try and find the fastest way possible of doing this, I've been pointed to a book called "Introduction to Algorithms", by T.H. cormen, C.E,Leisserson, R.L. Rivest.

    And that's given by a psuedocode algorithm that is a factor of 3 faster than methods such as using a matrix inverse. The problem comes when I try to translate the code to C.

    I realise this is a pretty complex topic, and most people probably won;t even want to bother wading through the code below, but I'm sorta out of options . If anyone's got some suggestions on other approaches or good sources of relevant information, I'd REALLY love to hear about it.

    Anyway, this is the psuedo-code I have:
    Code:
    LUP-solve (L,U,Pi,b)
    n = rows[L]
    for i = 1 to n
      do y[i] = b[Pi[i]] - sum(j=1 to i-1) of L[i][j]y[i]
    for i = n downto 1
      do x[i] = (y[i]-sum(j=i+1 to n) of u[i][j]x[j])/u[i][i]
    return x
    
    LUP-Decomposition(A)
    n = rows[A]
    for i = 1 to n
      do Pi[i] = i
    for k = 1 to n-1
      do p = 0
      for i = k to n
        do if |A[i][k]| > p
          then p=|A[i][k]|
            k' = i
      if p=0
        then error "singular matrix"
      exchange Pi[k] <> Pi[k']
      for i=1 to n
        do exchange A[k][i] <> A[k'][i]
      for i=k+1 to n
        do A[i][k] = A[i][k]/A[k][k]
        for j =k+1 to n
          do A[i][j] = A[i][j] - A[i][k]A[k][j]
    if anyone actually has the book, it's a fair bit easier to understand there, with the characters and formatting I just can't get in a forum post.

    My attempt at this code is:

    Code:
    int LUPsolve(int rows,float A[rows][rows],float b[rows],float X[rows])
    {
      int P[rows];
      int n,m,p,q,i,j;
      int k=0,kd=0;
      int T=0;
      float t=0;
      float Y[rows];
      //decompose
      for(i=1;i<rows;i++)
        {
        P[i]=i;
        }
      for(k=1;k<rows-1;k++)
        {
        p=0;
        for(i=k;i<rows;i++)
          {
          t=A[i][k];
          if(A[i][k]<0)t=-1*t;
          if(t>p)
            {
            p=t;
            kd=i;
            }
          }
        if(p==0)
          {
          printf("singular Matrix, no possible solutions\n");
          return 1;
          }
        
        T=P[kd];
        P[kd]=P[k];
        P[k]=T;
        for(i=1;i<rows;i++)
          {
          t=A[kd][i];
          A[kd][i]=A[k][i];
          A[k][i]=t;
          }
        for(i=k+1;i<rows;i++)
          {
          A[i][k]=A[i][k]/A[k][k];
          for(j=k+1;j<rows;j++)
            {
            A[i][j]=A[i][j]-A[i][k]*A[k][j];
            }
          }
        }
        
      //now solve
      for(n=1;n<rows;n++)
        {
        t=0;
        for(m=1;m<n-1;m++)
          {
          t+=A[n][m]*Y[m];
          }
        Y[n]=b[P[n]]-t;
        }
      for(n=rows-1;n>=1;n--)
        {
        t=0;
        for(m=n+1;m<rows;m++)
          {
          t+=A[n][m]*X[m];
          }
        X[n]=(Y[n]-t)/A[n][n];
        }
      //X now contains the solution;
      return 0;
    }
    I've compressed the two functions into one here, though that's not a problem, in the initial version, the decomposition runs, then solve, I've done the same thing, just in one function call.

    Also, though the psuedo-code for LUP-solve doesn't show it (decom does) the matrices L and U are both combined into A, so that is not the problem with the code. The indexing from 1 rather than 0 may be a problem, but when I go through and edit it to index from 0, it doesn;t work either, so I just keep it like this to maintain continuity with the book.

    the input I'm tesing it with is:

    a={{0,0,0},{0,4,2},{0,3,-2}}
    b={0,6,8}
    rows=3

    done on paper, the output should be 2 and -11, but I actually get 2.64 & -2.28

    anyone got an idea what's wrong/any pointers to more resources/pre-coded functions/anything remotely useful. thisis something I really want to be able to do, but my own brain's pretty much worn out on this

    big thanks to anyone who's still reading by now
    -mark
    "never argue with an idiot, they will drag you down to their level and beat you with experience"

  2. #2
    Just Lurking Dave_Sinkula's Avatar
    Join Date
    Oct 2002
    Posts
    5,005
    I don't know if this will be of any use, but here is some C++ code that might be helpful to look at.
    Code:
    int LUPsolve(int rows,float A[rows][rows],float b[rows],float X[rows])
    {
      int P[rows];
    This compiles? Is this C++? C99?
    7. It is easier to write an incorrect program than understand a correct one.
    40. There are two ways to write error-free programs; only the third one works.*

  3. #3
    Registered User Markallen85's Avatar
    Join Date
    Nov 2002
    Posts
    53
    it compiles and runs in dev_c++ in C mode.

    I know using "rows" all over th plpace there is not really a great way of doing things, but if I leave them out I get compile errors, something to do with variable arrays.

    Big thanks for that link, lines 85-128 seem exactly what I'm looking for, shame I don;tunderstand all the differences between C/C++ but I should be able to make out enough of it to dig around for the problems. I've bookmarked that and I'll deal with it in the morenig when I'm a little more awake (11:20 Pm is not the best time for hits kinda stuff...).

    thanks again

    -mark
    "never argue with an idiot, they will drag you down to their level and beat you with experience"

  4. #4
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    Mark, can you explain what U is. In the pseudo-code:
    LUP-solve (L,U,Pi,b)

    Now it appears that LUP-decomposition is done first, then LUP-solve. But in the pseudocode, U is never calculated. I'm assuming it's not the b array
    b=[0,6,8]
    as this seems to be passed last.

    I know you said L and U are combined into A, but you're using A as if it's L in the pseudocode.

  5. #5
    Registered User Markallen85's Avatar
    Join Date
    Nov 2002
    Posts
    53
    the pseudocode comes with 11 pages of further information detailing what's actually going on, a bit too much to repeat here, but I'll try and summarize what I can, though I should warn that I'm terrible at explanations:

    L and U are triangular matrices, (sort of only half a matrix) for example:
    L=
    1 0 0
    0.6 1 0
    0.2 0.5 1

    U=
    5 6 3
    0 1.4 2.2
    0 0 -1.8

    notice that they both have lots of zeros in the top right/bottom left corner.

    the LUP decomposition returns the matrix A, which contains both of these matrices:

    5 6 3
    0.6 1.4 2.2
    0.2 0.5 -1.8

    both L and U can be reconstructed from A, as the only cells that are not passed directly are ones.

    if you look at the sum statements in the LUP-Solve function, each sum is only summing one part of A.

    for i = 1 to n
    sum(j=1 to i-1)

    implies that (for n=3) we sum cells

    i=1: no cells
    i=2: j=1
    i=3: j=1, j=2

    and for the second:

    i=3:no cells
    i=2:j=3
    i=1:j=2,j=3

    so if we write the matri, showing which sum read each cell:

    _22
    1_2
    11_

    they're both reading triangles from the same matrix and I think I just saw a potential problem...the diagonal in the middle isna't being read at all, even though it's got values in. I think I'll take a closer look at that, thanks for making me step through it in detail again

    -mark
    "never argue with an idiot, they will drag you down to their level and beat you with experience"

  6. #6
    Registered User Markallen85's Avatar
    Join Date
    Nov 2002
    Posts
    53
    I'm really getting sick of silly errors, managed to have one in here that's been there for weeks.... at least I've got it now, the correct code is:

    Code:
    int LUPsolve(int rows,float A[rows][rows],float b[rows],float X[rows])
    {
      int P[rows];
      int n,m,p,q,i,j;
      int k=0,kd=0;
      int T=0;
      float t=0;
      float Y[rows];
      //initialise X,Y
      for (n=1;n<rows;n++)
        {
        X[n]=0;
        Y[n]=0;
        }
      //decompose
      for(i=1;i<rows;i++)
        {
        P[i]=i;
        }
      for(k=1;k<rows-1;k++)
        {
        p=0;
        for(i=k;i<rows;i++)
          {
          t=A[i][k];
          if(t<0)t=-1*t;
          if(t>p)
            {
            p=t;
            kd=i;
            }
          }
        if(p==0)
          {
          printf("singular Matrix, no possible solutions\n");
          return 1;
          }
        
        T=P[kd];
        P[kd]=P[k];
        P[k]=T;
        for(i=1;i<rows;i++)
          {
          t=A[kd][i];
          A[kd][i]=A[k][i];
          A[k][i]=t;
          }
        for(i=k+1;i<rows;i++)
          {
          A[i][k]=A[i][k]/A[k][k];
          for(j=k+1;j<rows;j++)
            {
            A[i][j]=A[i][j]-A[i][k]*A[k][j];
            }
          }
        }
      //now solve
      for(n=1;n<rows;n++)
        {
        t=0;
        for(m=1;m<=n-1;m++)
          {
          t+=A[n][m]*Y[m];
          }
        Y[n]=b[P[n]]-t;
        }
      for(n=rows-1;n>=1;n--)
        {
        t=0;
        for(m=n+1;m<rows;m++)
          {
          t+=A[n][m]*X[m];
          }
        X[n]=(Y[n]-t)/A[n][n];
        }
      //X now contains the solution;
      return 0;
    }
    the error was in the third last for statement, changing

    for(m=1;m<n-1;m++)

    to

    for(m=1;m<=n-1;m++)

    now works fine, dunno how I missed that in the who-knows-how many times I worked through the code & psuedocode... silly errors seems to be a really big thing for me

    -mark
    "never argue with an idiot, they will drag you down to their level and beat you with experience"

  7. #7
    Registered User
    Join Date
    Oct 2001
    Posts
    2,934
    Cool program. I ran your code with the values above and got 2 and -1. That was a tough bug to spot.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Read coefficients of linear equations from file into 2d array
    By omaralqady in forum C++ Programming
    Replies: 6
    Last Post: 06-20-2009, 07:39 AM
  2. Simultaneous linear Equations
    By dvd4alll in forum C# Programming
    Replies: 8
    Last Post: 02-08-2008, 07:11 PM
  3. Manipulation and analysis of Mathematical Equations
    By WebmasterMattD in forum C++ Programming
    Replies: 13
    Last Post: 06-07-2003, 10:08 PM
  4. simultaneous equations
    By nic_elmo_13 in forum C Programming
    Replies: 3
    Last Post: 04-14-2003, 08:11 AM
  5. equations for a program
    By anthonye in forum C Programming
    Replies: 4
    Last Post: 06-19-2002, 04:38 AM