Jacobi iteration on C#?

Question

Dear Expert, Do someone have about the Jacobi iteration source code on C#? please dear, i really need them for my research project. TQ..;)

nur saadah · Answer

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

Alan · Answer

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 { for (iq=ip+1; iq } if (sm == 0.0) break; if (i else tresh = 0.0; for (ip=0; 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 } 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 } 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 bool result = __jacobi(a, w, v); if( result ) { for(i=0;i for(i=0;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 bool result = __jacobi(a, w, v); if( result ) { for(i=0;i for(i=0;i } return result; } } // end Jacobi class

nur saadah · Answer

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 { for (iq=ip+1; iq } if (sm == 0.0) break; if (i else tresh = 0.0; for (ip=0; 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 } 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 } 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 bool result = __jacobi(a, w, v); if( result ) { for(i=0;i for(i=0;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 bool result = __jacobi(a, w, v); if( result ) { for(i=0;i for(i=0;i } return result; }

Alan · Answer

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 Console.WriteLine("
Invalid number, please re-enter
"); } Console.WriteLine("
Enter variable co-efficients, separated by spaces :
"); a = new double[n,n]; for (int i = 0; i { do { Console.Write(" Equation {0} : ", i + 1); temp = Console.ReadLine(); tempArray = temp.Split(new char[]{' '}, StringSplitOptions.RemoveEmptyEntries); if (tempArray.Length != n) { Console.WriteLine("
Invalid number of coefficients, please re-enter
"); isValidNumber = false; continue; } for (int j = 0; j { isValidNumber = double.TryParse(tempArray[j], out value); if (!isValidNumber) { Console.WriteLine("
Line contains an invalid number, please re-enter whole line
"); break; } else if (j == i && value == 0) { Console.WriteLine("
Main diagonal cannot contain zero co-efficients, please re-enter whole line
"); isValidNumber = false; break; } else { a[i,j] = value; } } } while(!isValidNumber); } Console.WriteLine("
Enter constant values, separated by spaces
"); 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("
Invalid number of constant values, please re-enter
"); isValidNumber = false; continue; } for (int i = 0; i { isValidNumber = double.TryParse(tempArray[i], out value); if (!isValidNumber) { Console.WriteLine("
Line contains an invalid number, please re-enter whole line
"); break; } else { b[i] = value; } } } while(!isValidNumber); Console.WriteLine("
Enter initial approximations, separated by spaces
"); 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("
Invalid number of approximations, please re-enter
"); isValidNumber = false; continue; } for (int i = 0; i { isValidNumber = double.TryParse(tempArray[i], out value); if (!isValidNumber) { Console.WriteLine("
Line contains an invalid number, please re-enter whole line
"); 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("
Invalid number, please re-enter
"); } 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 Console.WriteLine("
Invalid number, please re-enter
"); } Console.WriteLine(); iterations = max; withinTol = false; for (int iteration = 1 ; iteration { finished = true; for (int i = 0; i { x[i] = b[i]; for(int j = 0; 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 :
"); for (int i = 1; i { Console.Write(" x{0} = {1:F5}", i, x[i - 1]); // display to 5 dp Console.WriteLine(); } Console.WriteLine("
Number of iterations : {0}", iterations); Console.WriteLine("Approximations are within tolerance : {0}", withinTol); Console.Write("
Press any key to exit program
"); Console.ReadKey(); } }

Alan · Answer

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/

Jacobi iteration on C#?

Answers (5)