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[ipip] = 1.0; } for (ip=0; ip { b[ip] = w[ip] = a[ipip]; 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[ipiq]); if (i > 4 && (Math.Abs(w[ip])+g) == Math.Abs(w[ip]) && (Math.Abs(w[iq])+g) == Math.Abs(w[iq])) { a[ipiq] = 0.0; } else if (Math.Abs(a[ipiq]) > tresh) { h = w[iq] - w[ip]; if ( (Math.Abs(h)+g) == Math.Abs(h)) t = (a[ipiq]) / h; else { theta = 0.5*h / (a[ipiq]); 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[ipiq]; z[ip] -= h; z[iq] += h; w[ip] -= h; w[iq] += h; a[ipiq]=0.0; for (j=0;j { g=a[jip]; h=a[jiq]; a[jip]=g-s*(h+g*tau); a[jiq]=h+s*(g-h*tau); } for (j=ip+1;j { g=a[ipj]; h=a[jiq]; a[ipj]=g-s*(h+g*tau); a[jiq]=h+s*(g-h*tau); } for (j=iq+1; j { g=a[ipj]; h=a[iqj]; a[ipj]=g-s*(h+g*tau); a[iqj]=h+s*(g-h*tau); } for (j=0; j { g=v[jip]; h=v[jiq]; v[jip]=g-s*(h+g*tau); v[jiq]=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[ij]; v[ij] = v[ik]; v[ik] = tmp; } } } // VTK addition to original Numerical Recipes code: // insure eigenvector consistency (i.e. Jacobi can compute // vectors that are negative of one another (.707.7070) and // (-.707-.7070). 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[33]; double[] v = new double[33]; double[] w = new double[3]; int ij; 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[44]; double[] v = new double[44]; double[] w = new double[4]; int ij; 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.cxxv 1.2 2001/09/24 21:57:17 garland Exp $ ************************************************************************/ #include gfx\gfx.h #include gfx\mat3.h #include gfx\mat4.h #define ROT(aijkl) 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(ajipjiq); } for (j=ip+1;j { ROT(aipjjiq); } for (j=iq+1; j { ROT(aipjiqj); } for (j=0; j { ROT(vjipjiq); } } } } 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.7070) and // (-.707-.7070). 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 ij; 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 ij; 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[nn]; 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[ij] = 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[ii]; 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/