Jacobi iteration on C#?

Tq so much..Alan..
I get your idea..

;)

Here's my attempt at a C# translation of that C++ code. I've changed it as little as possible so as not to disturb the algorithms but notice in particular that:

1. As C# doesn't have support for macros, I've simply expanded ROT and made MAX_ROTATIONS into a private constant.

2. As the .NET framework doesn't have support for either matrices or vectors, I've represented the former by 2 dimensional and the latter by single dimensional double arrays.

3. C# doesn't, of course, have templates and there was no point in declaring __jacobi as a generic method as it simply won't compile.

4. As there's no way to declare a parameter which takes a fixed size array in C#, I've had to rename the two 'eigen' functions eigen3 and eigen4 otherwise the compiler would have had no way of distinguishing their signatures. I've left the eig_vecs parameter as an array of single dimensional arrays rather than a 'normal' 2 dimensional array in C#.

5. As C# doesn't support 'standalone' functions, I've made all three methods public static and put them inside a public static class called Jacobi. 

public static class Jacobi
{

private const int MAX_ROTATIONS = 60;

public static bool __jacobi(double[,] a, double[] w, double[,] v)
{
    int N = w.Length;
    int i, j, k, iq, ip;
    double tresh, theta, tau, t, sm, s, h, g, c;
    double tmp;
    double[] b = new double[N];
    double[] z = new double[N];

    // initialize
    for (ip=0; ip    {
        for (iq=0; iq        v[ip,ip] = 1.0;
    }
    for (ip=0; ip    {
        b[ip] = w[ip] = a[ip,ip];
        z[ip] = 0.0;
    }

    // begin rotation sequence
    for (i=0; i    {
        sm = 0.0;
        for (ip=0; ip<2; ip++)
        {
            for (iq=ip+1; iq        }
        if (sm == 0.0) break;

        if (i < 4) tresh = 0.2*sm/(9);
        else tresh = 0.0;

        for (ip=0; ip<2; ip++)
        {
            for (iq=ip+1; iq            {
                g = 100.0 * Math.Abs(a[ip,iq]);
                if (i > 4 && (Math.Abs(w[ip])+g) == Math.Abs(w[ip])
                    && (Math.Abs(w[iq])+g) == Math.Abs(w[iq]))
                {
                    a[ip,iq] = 0.0;
                }
                else if (Math.Abs(a[ip,iq]) > tresh)
                {
                    h = w[iq] - w[ip];
                    if ( (Math.Abs(h)+g) == Math.Abs(h)) t = (a[ip,iq]) / h;
                    else
                    {
                        theta = 0.5*h / (a[ip,iq]);
                        t = 1.0 / (Math.Abs(theta) + Math.Sqrt(1.0+theta*theta));
                        if (theta < 0.0) t = -t;
                    }
                    c = 1.0 / Math.Sqrt(1+t*t);
                    s = t*c;
                    tau = s/(1.0+c);
                    h = t*a[ip,iq];
                    z[ip] -= h;
                    z[iq] += h;
                    w[ip] -= h;
                    w[iq] += h;
                    a[ip,iq]=0.0;
                    for (j=0;j                    {
                        g=a[j,ip];
                        h=a[j,iq];
                        a[j,ip]=g-s*(h+g*tau);
                        a[j,iq]=h+s*(g-h*tau);
                    }
                    for (j=ip+1;j                    {
                        g=a[ip,j];
                        h=a[j,iq];
                        a[ip,j]=g-s*(h+g*tau);
                        a[j,iq]=h+s*(g-h*tau);
                    }
                    for (j=iq+1; j                    {
                        g=a[ip,j];
                        h=a[iq,j];
                        a[ip,j]=g-s*(h+g*tau);
                        a[iq,j]=h+s*(g-h*tau);
                    }
                    for (j=0; j                    {
                        g=v[j,ip];
                        h=v[j,iq];
                        v[j,ip]=g-s*(h+g*tau);
                        v[j,iq]=h+s*(g-h*tau);
                    }
                }
            }
        }

        for (ip=0; ip        {
            b[ip] += z[ip];
            w[ip] = b[ip];
            z[ip] = 0.0;
        }
    }

    if ( i >= MAX_ROTATIONS )
        return false;

    // sort eigenfunctions
    for (j=0; j    {
        k = j;
        tmp = w[k];
        for (i=j; i        {
            if (w[i] >= tmp)
            {
                k = i;
                tmp = w[k];
            }
        }
        if (k != j)
        {
            w[k] = w[j];
            w[j] = tmp;
            for (i=0; i            {
                tmp = v[i,j];
                v[i,j] = v[i,k];
                v[i,k] = tmp;
            }
        }
    }

    // VTK addition to original Numerical Recipes code:
    //    insure eigenvector consistency (i.e., Jacobi can compute
    //    vectors that are negative of one another (.707,.707,0) and
    //    (-.707,-.707,0). This can reek havoc in
    //    hyperstreamline/other stuff. We will select the most
    //    positive eigenvector.
    int numPos;
    for (j=0; j    {
        for (numPos=0, i=0; i= 0.0 ) numPos++;
        if ( numPos < 2 ) for(i=0; i    }

    return true;
}

public static bool eigen3(double[,] m, double[] eig_vals, double[][] eig_vecs)
{
    double[,] a = new double[3,3];
    double[,] v = new double[3,3];
    double[] w = new double[3];
    int i,j;

    for(i=0;i<3;i++) for(j=0;j<3;j++) a[i,j] = m[i,j];

    bool result = __jacobi(a, w, v);
    if( result )
    {
        for(i=0;i<3;i++) eig_vals[i] = w[i];

        for(i=0;i<3;i++) for(j=0;j<3;j++) eig_vecs[i][j] = v[j,i];
    }

    return result;
}

public static bool eigen4(double[,] m, double[] eig_vals, double[][] eig_vecs)
{
    double[,] a = new double[4,4];
    double[,] v = new double[4,4];
    double[] w = new double[4];
    int i,j;

    for(i=0;i<4;i++) for(j=0;j<4;j++) a[i,j] = m[i,j];

    bool result = __jacobi(a, w, v);
    if( result )
    {
        for(i=0;i<4;i++) eig_vals[i] = w[i];

        for(i=0;i<4;i++) for(j=0;j<4;j++) eig_vecs[i][j] = v[j,i];
    }

    return result;
}

} // end Jacobi class

thanks for your response.;)

for your information, i have jacobi iteration source code in C language.
The problem is, i face problem to convert it to C# because some of their function cannot
use in C# such as  preprocessor macros and so on.until now, i still try to solve it but still have problem.i need ur help.. TQ;)

the source code is below:

/************************************************************************

Jacobi iteration routines for computing eigenvalues/eigenvectors.

NOTE: This code was adapted from VTK source code (vtkMath.cxx) which
seems to have been adapted directly from Numerical Recipes in C.

$Id: jacobi.cxx,v 1.2 2001/09/24 21:57:17 garland Exp $

************************************************************************/
 
#include "gfx\gfx.h"
#include "gfx\mat3.h"
#include "gfx\mat4.h"

#define ROT(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau)

#define MAX_ROTATIONS 60

// Description:
// Jacobi iteration for the solution of eigenvectors/eigenvalues of a NxN
// real symmetric matrix. Square NxN matrix a; output eigenvalues in w;
// and output eigenvectors in v. Resulting eigenvalues/vectors are sorted
// in decreasing order; eigenvectors are normalized.
 
template
bool __jacobi(Mat &a, Vec &w, Mat &v)
{
    const int N = Vec::dim();

    int i, j, k, iq, ip;
    double tresh, theta, tau, t, sm, s, h, g, c;
    double tmp;
    Vec b, z;

    // initialize
    for (ip=0; ip     {
        for (iq=0; iq         v[ip][ip] = 1.0;
    }
    for (ip=0; ip     {
        b[ip] = w[ip] = a[ip][ip];
        z[ip] = 0.0;
    }

    // begin rotation sequence
    for (i=0; i     {
        sm = 0.0;
        for (ip=0; ip<2; ip++)
        {
            for (iq=ip+1; iq         }
        if (sm == 0.0) break;

        if (i < 4) tresh = 0.2*sm/(9);
        else tresh = 0.0;

        for (ip=0; ip<2; ip++)
        {
            for (iq=ip+1; iq             {
                g = 100.0*fabs(a[ip][iq]);
                if (i > 4 && (fabs(w[ip])+g) == fabs(w[ip])
                    && (fabs(w[iq])+g) == fabs(w[iq]))
                {
                    a[ip][iq] = 0.0;
                }
                else if (fabs(a[ip][iq]) > tresh)
                {
                    h = w[iq] - w[ip];
                    if ( (fabs(h)+g) == fabs(h)) t = (a[ip][iq]) / h;
                    else
                    {
                        theta = 0.5*h / (a[ip][iq]);
                        t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
                        if (theta < 0.0) t = -t;
                    }
                    c = 1.0 / sqrt(1+t*t);
                    s = t*c;
                    tau = s/(1.0+c);
                    h = t*a[ip][iq];
                    z[ip] -= h;
                    z[iq] += h;
                    w[ip] -= h;
                    w[iq] += h;
                    a[ip][iq]=0.0;
                    for (j=0;j                     {
                        ROT(a,j,ip,j,iq);
                    }
                    for (j=ip+1;j                     {
                        ROT(a,ip,j,j,iq);
                    }
                    for (j=iq+1; j                     {
                        ROT(a,ip,j,iq,j);
                    }
                    for (j=0; j                     {
                        ROT(v,j,ip,j,iq);
                    }
                }
            }
        }

        for (ip=0; ip         {
            b[ip] += z[ip];
            w[ip] = b[ip];
            z[ip] = 0.0;
        }
    }

    if ( i >= MAX_ROTATIONS )
        return false;

    // sort eigenfunctions
    for (j=0; j     {
        k = j;
        tmp = w[k];
        for (i=j; i         {
            if (w[i] >= tmp)
            {
                k = i;
                tmp = w[k];
            }
        }
        if (k != j)
        {
            w[k] = w[j];
            w[j] = tmp;
            for (i=0; i             {
                tmp = v[i][j];
                v[i][j] = v[i][k];
                v[i][k] = tmp;
            }
        }
    }

    // VTK addition to original Numerical Recipes code:
    //    insure eigenvector consistency (i.e., Jacobi can compute
    //    vectors that are negative of one another (.707,.707,0) and
    //    (-.707,-.707,0). This can reek havoc in
    //    hyperstreamline/other stuff. We will select the most
    //    positive eigenvector.
    int numPos;
    for (j=0; j     {
        for (numPos=0, i=0; i= 0.0 ) numPos++;
        if ( numPos < 2 ) for(i=0; i     }

    return true;
}

#undef ROT
#undef MAX_ROTATIONS

bool eigen(const Mat3& m, Vec3& eig_vals, Vec3 eig_vecs[3])
{
    Mat3 a, v;  Vec3 w;
    int i,j;

    for(i=0;i<3;i++) for(j=0;j<3;j++) a[i][j] = m(i,j);

    bool result = __jacobi(a, w, v);
    if( result )
    {
        for(i=0;i<3;i++) eig_vals[i] = w[i];

        for(i=0;i<3;i++) for(j=0;j<3;j++) eig_vecs[i][j] = v[j][i];
    }

    return result;
}

bool eigen(const Mat4& m, Vec4& eig_vals, Vec4 eig_vecs[4])
{
    Mat4 a, v;  Vec4 w;
    int i,j;

    for(i=0;i<4;i++) for(j=0;j<4;j++) a[i][j] = m(i,j);

    bool result = __jacobi(a, w, v);
    if( result )
    {
        for(i=0;i<4;i++) eig_vals[i] = w[i];

        for(i=0;i<4;i++) for(j=0;j<4;j++) eig_vecs[i][j] = v[j][i];
    }

    return result;
}

As the algorithm is straightforward, I thought I'd have a go at implementing  the Jacobi iteration method in C#. The console application below (.NET 2.0 or later) seems to be working OK but, if you want to use it, I'd test it thoroughly first:

using System;

class Jacobi
{
    static void Main()
    {
       Console.Clear();    
       int n ;      // number of equations/variables
       double[,] a; // co-efficients of variables (on LHS)
       double[] b;  // constant values (on RHS)
       double[] x0; // previous approximation to variable values
       double[] x;  // current approximation to variable values
       double[] diff; // absolute difference between approximations
       double tol;  // tolerance
       int max;     // maximum number of iterations      
       int iterations; // actual number of iterations required
       bool withinTol; // whether the results are within the tolerance
      
       bool isValidNumber;
       string temp;
       string[] tempArray;
       double value;
       bool finished;

       while(true)
       {
           Console.Write("Enter number of equations/variables (2 to 20): ");
           isValidNumber = int.TryParse(Console.ReadLine(), out n);
           if(isValidNumber && n > 1 && n < 21) break;
           Console.WriteLine("\nInvalid number, please re-enter\n");
       } 

       Console.WriteLine("\nEnter variable co-efficients, separated by spaces :\n");
       a = new double[n,n];
   
      
       for (int i = 0; i < n ; i++)
       {
           do
           {
               Console.Write(" Equation {0} : ", i + 1);
               temp = Console.ReadLine();
               tempArray = temp.Split(new char[]{' '}, StringSplitOptions.RemoveEmptyEntries);
               if (tempArray.Length != n)
               {
                   Console.WriteLine("\nInvalid number of coefficients, please re-enter\n");
                   isValidNumber = false;
                   continue;
               }
               for (int j = 0; j < n; j++)
               {
                   isValidNumber = double.TryParse(tempArray[j], out value);
                   if (!isValidNumber)
                   {
                       Console.WriteLine("\nLine contains an invalid number, please re-enter whole line\n");
                       break;
                   }
                   else if (j == i && value == 0)
                   {
                       Console.WriteLine("\nMain diagonal cannot contain zero co-efficients, please re-enter whole line\n");
                       isValidNumber = false;
                       break;
                   }
                   else
                   {
                       a[i,j] = value;
                   }
               }
           }
           while(!isValidNumber);
       }
       
       Console.WriteLine("\nEnter constant values, separated by spaces\n");
       b = new double[n];         
      
       do
       {
           Console.Write(" For all equations : ");
           temp = Console.ReadLine();
           tempArray = temp.Split(new char[]{' '}, StringSplitOptions.RemoveEmptyEntries);
           if (tempArray.Length != n)
           {
               Console.WriteLine("\nInvalid number of constant values, please re-enter\n");
               isValidNumber = false;
               continue;
           }
           for (int i = 0; i < n; i++)
           {
               isValidNumber = double.TryParse(tempArray[i], out value);
               if (!isValidNumber)
               {
                   Console.WriteLine("\nLine contains an invalid number, please re-enter whole line\n");
                   break;
               }
               else
               {
                   b[i] = value;
               }
           }         
 
       }
       while(!isValidNumber);
 
       Console.WriteLine("\nEnter initial approximations, separated by spaces\n");
       x0 = new double[n];
       x = new double[n];
       diff = new double[n];

       do
       {
           Console.Write(" For all variables : ");
           temp = Console.ReadLine();
           tempArray = temp.Split(new char[]{' '}, StringSplitOptions.RemoveEmptyEntries);
           if (tempArray.Length != n)
           {
               Console.WriteLine("\nInvalid number of approximations, please re-enter\n");
               isValidNumber = false;
               continue;
           }
           for (int i = 0; i < n; i++)
           {
               isValidNumber = double.TryParse(tempArray[i], out value);
               if (!isValidNumber)
               {
                   Console.WriteLine("\nLine contains an invalid number, please re-enter whole line\n");
                   break;
               }
               else
               {
                   x0[i] = value;
               }
           }          
           
       }
       while(!isValidNumber);
 
       Console.WriteLine();

       while(true)
       {
           Console.Write("Enter tolerance ( > 0) for all variables : ");
           isValidNumber = double.TryParse(Console.ReadLine(), out tol);
           if(isValidNumber && tol > 0) break;
           Console.WriteLine("\nInvalid number, please re-enter\n");
       }

       Console.WriteLine();   

       while(true)
       {
           Console.Write("Enter maximum number of iterations (5 to 99) : ");
           isValidNumber = int.TryParse(Console.ReadLine(), out max);
           if(isValidNumber &&  max > 4 && max < 100) break;
           Console.WriteLine("\nInvalid number, please re-enter\n");
       }

       Console.WriteLine();
       iterations  = max;
       withinTol  = false;
          
       for (int iteration = 1 ; iteration <= max; iteration++)
       {
           finished = true;
           for (int i = 0; i < n; i++)
           {
               x[i] = b[i];
               for(int j = 0; j < n; j++)            
               {
                   if (j == i) continue;
                   x[i] -= a[i, j] * x0[j];                  
               }            
               x[i] /= a[i,i];
               diff[i] = Math.Abs(x[i] - x0[i]);
               if (finished && diff[i] > tol) finished = false;
           }
           if (finished)
           {
               iterations = iteration;
               withinTol = true;
               break;
           }         
           Array.Copy(x, x0, n);
       }
                  
       Console.WriteLine("The approximate values of the variables are :\n");
       for (int i = 1; i <= n; i++)
       {
           Console.Write(" x{0} = {1:F5}", i, x[i - 1]); // display to 5 dp
           Console.WriteLine(); 
       }
 
       Console.WriteLine("\nNumber of iterations : {0}", iterations);
       Console.WriteLine("Approximations are within tolerance : {0}", withinTol);
     
       Console.Write("\nPress any key to exit program\n");
       Console.ReadKey();
       
   }
}

I can't find any C# code at all which deals with the Jacobi iteration method.

However, it might be worth having a look at the free linear algebra library mapack for .net which does, I believe, include an equation solver:

http://www.aisto.com/roeder/dotnet/