fork download
copy 
  1. #include <stdio.h>
  2. #include <math.h>
  3. #define m1 201
  4. #define n1 101
  5. #define Max(a,b) (a>b ? a:b)
  6. int HQRdecomp(int m, int n, double eps, double a[][n1], double d[]){
  7. 	double r,s,p;
  8. 	int i,j,k;
  9. 	for(k=1;k<=n;k++){
  10. 		s=0;
  11. 		for(i=k;i<=m;i++){
  12. 			s+=a[i][k]*a[i][k];
  13. 		}
  14. 		s=sqrt(s);
  15. 		if(s<eps){
  16. 			return k;
  17. 		}
  18. 		d[k]=(a[k][k]>0)? -s:s;
  19. 		r=1/sqrt(d[k]*(d[k]-a[k][k]));
  20. 		a[k][k]=(a[k][k]-d[k])*r;
  21. 		for(i=k+1;i<=m;i++){
  22. 			a[i][k]*=r;
  23. 		}
  24. 		for(j=k+1;j<=n;j++){
  25. 			p=0;
  26. 			for(i=k;i<=m;i++){
  27. 				p+=a[i][j]*a[i][k];
  28. 			}
  29. 			for(i=k;i<=m;i++){
  30. 				a[i][j]-=p*a[i][k];
  31. 			}
  32. 		}
  33. 	}
  34. 	return 0;
  35. }
  36. int HQRsolve(double a[][n1],double d[],int m,int n,double b[]){
  37. 	double p;
  38. 	int i,j,k;
  39. 	for(k=1;k<=n;k++){
  40. 		p=0;
  41. 		for(i=k;i<=m;i++){
  42. 			p+=b[i]*a[i][k];
  43. 		}
  44. 		for(i=k;i<=m;i++){
  45. 			b[i]-=p*a[i][k];
  46. 		}
  47. 	}
  48. 	for(i=n;i>=1;i--){
  49. 		for(j=i+1;j<=n;j++){
  50. 			b[i]-=a[i][j]*b[j];
  51. 		}
  52. 		b[i]=b[i]/d[i];
  53. 	}
  54. 	return 0;
  55. }
  56. int main(void) {
  57. 	double a[m1][n1],d[m1],b[m1],eps;
  58. 	int i,j,m,n,cnd;
  59. 	m=4;
  60. 	n=2;
  61. 	a[1][1]=1;
  62. 	a[1][2]=1;
  63. 	a[2][1]=1;
  64. 	a[2][2]=3;
  65. 	a[3][1]=1;
  66. 	a[3][2]=6;
  67. 	a[4][1]=1;
  68. 	a[4][2]=9;
  69. 	b[1]=3;
  70. 	b[2]=5;
  71. 	b[3]=7;
  72. 	b[4]=12;
  73.  
  74. 	eps=1.0e-14;
  75. 	cnd=HQRdecomp(m,n,eps,a,d);
  76. 	printf("cnd=%d\n",cnd);
  77. 	HQRsolve(a,d,m,n,b);
  78. 	for(i=1;i<=n;i++){
  79. 		printf("b[%2d]=%10.7f\n",i,b[i]);
  80. 	}
  81. 	return 0;
  82. }
  83.  
Success #stdin #stdout 0s 5692KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
cnd=0
b[ 1]= 1.6122449
b[ 2]= 1.0816327
https://ideone.com/mJybPA
language:
C (gcc 8.3)
created:
2 years ago
visibility:
public

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!