fork download
copy 
  1. using System;
  2. using System.Globalization;
  3. namespace NamespaceName
  4. {
  5. 	public class ClassName
  6. 	{
  7. 		const int maxiter = 1000;
  8. 		const double eps = 1e-12;
  9. 		public static void PrintMatrix(int m,int n,double[,] A)
  10. 		{
  11. 			NumberFormatInfo nfi = new NumberFormatInfo();
  12. 			nfi.NumberDecimalDigits = 12;
  13. 			double p = 0;
  14. 			for(int i=0;i<m;i++)
  15. 			{
  16. 				for(int j=0;j < n; j++)
  17. 					if(Math.Abs(A[i,j]) > eps)
  18. 						Console.Write("{0} , ",A[i,j].ToString("N",nfi));
  19. 					else
  20. 						Console.Write("{0} , ",p.ToString("N",nfi)); 
  21. 				Console.WriteLine();
  22. 			}
  23. 			Console.WriteLine();
  24. 		}
  25. 		public static double Hypot(double a,double b)
  26. 		{
  27. 			double absa,absb;
  28. 			absa = Math.Abs(a);
  29. 			absb = Math.Abs(b);
  30. 			if(absa > absb) return absa * Math.Sqrt(1 + ((double) absb / absa) * ((double) absb / absa));
  31. 			else return (absb == 0 ? 0 : absb * Math.Sqrt(1 + ((double) absa / absb) * ((double) absa / absb)));
  32. 		}
  33. 		public static void ToHessenberg(int n,double[,] A)
  34. 		{
  35. 			for(int i = 0;i < n-2;i++)
  36. 			{
  37. 				int k = i + 1;
  38. 				double maxA = Math.Abs(A[k,i]);
  39. 				for(int j= i + 1; j < n; j++)
  40. 				   if(Math.Abs(A[j,i]) > maxA)
  41. 				   {
  42. 					   k  = j;
  43. 					   maxA = Math.Abs(A[k,i]);
  44. 				   }
  45. 				if(maxA > 0)
  46. 				{
  47. 					for(int j = 0; j < n; j++)
  48. 					{
  49. 						double temp = A[k,j];
  50. 						A[k,j] = A[i + 1,j];
  51. 						A[i + 1,j] = temp;
  52. 					}
  53. 					for(int j=0;j < n;j++)
  54. 					{
  55. 						double temp = A[j,k];
  56. 						A[j,k] = A[j, i + 1];
  57. 						A[j,i + 1] = temp;
  58. 					}
  59. 					for(int j = i + 2;j < n;j++)
  60. 					{
  61. 						double m = (double)A[j,i]/A[i + 1,i];
  62. 						for(int l = 0; l < n;l++)
  63. 							A[j,l] -= m * A[i + 1,l];
  64. 						for(int l = 0;l < n;l++)
  65. 							A[l,i + 1] += m * A[l,j];
  66. 					}
  67. 				}
  68. 			}
  69. 		}
  70. 		public static void QRdecomp(int m,int n,double[,] A,double[,] Q)
  71. 		{
  72. 			for(int i = 0; i < m; i++)
  73. 				for(int j = 0;j < m;j++)
  74. 					Q[i,j] = (i == j ? 1 : 0);
  75. 				int min = (m < n ? m : n);
  76. 				for(int i = 0; i < min; i++)
  77. 					for(int j = i + 1; j < m; j++)
  78. 						if(A[j,i] != 0)
  79. 						{
  80. 							double r = Hypot(A[i,i],A[j,i]);
  81. 							double c = (double) A[i,i] / r;
  82. 							double s = (double) A[j,i] / r;
  83. 							for(int k = 0; k < n; k++)
  84. 							{
  85. 								double temp = A[i,k];
  86. 								A[i,k] = c * A[i,k] + s * A[j,k];
  87. 								A[j,k] = -s * temp + c * A[j,k];
  88. 							}
  89. 							for(int k = 0;k < m;k++)
  90. 							{
  91. 								double temp = Q[k,i];
  92. 								Q[k,i] = c * Q[k,i] + s * Q[k,j];
  93. 								Q[k,j] = -s * temp + c * Q[k,j];
  94. 							}
  95. 						}
  96.  
  97. 		}
  98.  
  99. 		public static void MultiplyMatrix(int m,int n,int p,double[,] A,double[,] B,double[,] C)
  100. 		{
  101. 			for(int i = 0;i < m;i++)
  102. 				for(int k=0;k < p;k++)
  103. 				{
  104. 					C[i,k] = 0;
  105. 					for(int j=0;j < n;j++)
  106. 						C[i,k] += A[i,j] * B[j,k];
  107. 				}
  108. 		}
  109. 		public static void CopyMatrix(int m,int n,double[,] A,double[,] B)
  110. 		{
  111. 			for(int i=0;i<m;i++)
  112. 				for(int j=0;j<n;j++)
  113. 					B[i,j] = A[i,j];
  114. 		}
  115. 		public static void Main(string[] args)
  116. 		{
  117. 			int n,m;
  118. 			double[,] A;
  119. 			double[,] Q;
  120. 			double[,] R;
  121. 			Console.WriteLine("Obliczanie przyblizonych wartosci wlasnych metoda rozkladu QR");
  122. 				Console.WriteLine("Podaj wymiar macierzy ktorej wartosci wlasne chcesz policzyc");
  123. 				int.TryParse(Console.ReadLine(),out n);
  124. 				Console.WriteLine("Podaj gorny zakres przedzialu do ktorego beda nalezec elementy macierzy");
  125. 				int.TryParse(Console.ReadLine(),out m);
  126. 				A = new double[n,n];
  127. 				Q = new double[n,n];
  128. 				R = new double[n,n];
  129. 				Random rand = new Random();
  130. 				for(int i=0;i<n;i++)
  131. 					for(int j = 0;j < n;j++)
  132. 						A[i,j] = (1 - 2 * rand.Next(2)) * rand.Next(m);
  133. 				Console.WriteLine("Wylosowana macierz");
  134. 				PrintMatrix(n,n,A);
  135. 				ToHessenberg(n,A);
  136. 				Console.WriteLine("Macierz po sprowadzeniu do postaci Hessenberga");
  137. 				PrintMatrix(n,n,A);
  138. 				for(int i = 0;i<maxiter;i++)
  139. 				{
  140. 					QRdecomp(n,n,A,Q);
  141. 					CopyMatrix(n,n,A,R);
  142. 					MultiplyMatrix(n,n,n,R,Q,A);
  143. 				}
  144. 				Console.WriteLine("Macierz po iteracyjnym sprowadzeniu do postaci Schura");
  145. 				PrintMatrix(n,n,A);
  146. 				int k = 0;
  147. 				NumberFormatInfo nfi = new NumberFormatInfo();
  148. 				nfi.NumberDecimalDigits = 12;
  149. 				while(k < n)
  150. 				{
  151. 					if(k + 1 < n && Math.Abs(A[k+1,k]) > eps)
  152. 					{
  153. 						/* Tutaj zalozylem ze w postaci Schura klatka 2x2 pojawia sie
  154. 						 na diagonali gdy macierz posiada zespolone wartosci wlasne
  155. 						 */
  156. 						double p = 0.5 * (A[k,k] + A[k + 1,k + 1]);
  157. 						double q = A[k,k] * A[k+1,k+1] - A[k,k+1] * A[k+1,k];
  158. 						double d = Math.Sqrt(Math.Abs(q - p * p));
  159. 						Console.WriteLine("x[{0}] = {1} - {2}*i",k,p.ToString("N",nfi),d.ToString("N",nfi));
  160. 						Console.WriteLine("x[{0}] = {1} + {2}*i",k+1,p.ToString("N",nfi),d.ToString("N",nfi));
  161. 						k += 2;
  162. 					}
  163. 					else
  164. 					{
  165. 						Console.WriteLine("x[{0}] = {1}",k,A[k,k].ToString("N",nfi));
  166. 						k++;
  167. 					}
  168. 				}		
  169. 		}
  170. 	}
  171. }
  172.  
Success #stdin #stdout 0.04s 25284KB
comments (?)
stdin
copy 
6
100
stdout
copy 
Obliczanie przyblizonych wartosci wlasnych metoda rozkladu QR
Podaj wymiar macierzy ktorej wartosci wlasne chcesz policzyc
Podaj gorny zakres przedzialu do ktorego beda nalezec elementy macierzy
Wylosowana macierz
50.000000000000 , -4.000000000000 , 36.000000000000 , -62.000000000000 , 26.000000000000 , -55.000000000000 , 
-18.000000000000 , -19.000000000000 , 9.000000000000 , 87.000000000000 , 28.000000000000 , 62.000000000000 , 
60.000000000000 , 41.000000000000 , 94.000000000000 , 50.000000000000 , -38.000000000000 , -88.000000000000 , 
-70.000000000000 , -41.000000000000 , -52.000000000000 , 96.000000000000 , -62.000000000000 , 22.000000000000 , 
-8.000000000000 , 39.000000000000 , 30.000000000000 , -89.000000000000 , 91.000000000000 , -51.000000000000 , 
-37.000000000000 , -22.000000000000 , 0.000000000000 , -29.000000000000 , 73.000000000000 , -55.000000000000 , 

Macierz po sprowadzeniu do postaci Hessenberga
50.000000000000 , -119.985714285714 , -34.428814698633 , 45.988245403736 , -28.474520843633 , -55.000000000000 , 
-70.000000000000 , 134.571428571429 , -1.967376204347 , -30.866002055952 , -31.210191662547 , 22.000000000000 , 
0.000000000000 , -136.622448979592 , 12.606021574213 , 17.595722952219 , 19.872323355780 , -53.514285714286 , 
0.000000000000 , 0.000000000000 , -170.402570335754 , 61.290052786869 , -18.443749698397 , -86.558082007618 , 
0.000000000000 , 0.000000000000 , 0.000000000000 , 65.457768486961 , 16.755848740297 , -24.394519807795 , 
0.000000000000 , 0.000000000000 , 0.000000000000 , 0.000000000000 , -58.394467123537 , -18.223351672808 , 

Macierz po iteracyjnym sprowadzeniu do postaci Schura
135.121802693143 , 122.544269149596 , -21.144163199209 , 11.909137803369 , -14.569156253241 , -59.969193765852 , 
-31.158058040203 , 149.267345057168 , 21.542275595957 , -86.602047931168 , 152.994725309717 , -4.435405604538 , 
0.000000000000 , 0.000000000000 , -36.618093965312 , 70.318210019517 , -94.022615149702 , -6.686945621790 , 
0.000000000000 , 0.000000000000 , -46.866188359712 , -53.549371565518 , 53.760852706482 , 15.178537351671 , 
0.000000000000 , 0.000000000000 , 0.000000000000 , 16.962298883562 , 67.424836305601 , 29.199783169510 , 
0.000000000000 , 0.000000000000 , 0.000000000000 , 0.000000000000 , 0.000000000000 , -4.646518525077 , 

x[0] = 142.194573875155 - 61.385807467722*i
x[1] = 142.194573875155 + 61.385807467722*i
x[2] = -45.083732765415 - 56.779216581410*i
x[3] = -45.083732765415 + 56.779216581410*i
x[4] = 67.424836305601
x[5] = -4.646518525077
https://ideone.com/c7vgTD
language:
C# (gmcs 5.20.1)
created:
2 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!