CMRH.m 3.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136
  1. function [x,tnorm,flag]=CMRH(A,b,tol,maxiter,x0,iprint)
  2. %
  3. % [x,tnorm,flag] = (A,b,tol,maxiter,x0,iprint)
  4. %
  5. % This function solves the linear system Ax=b using the CMRH method
  6. %
  7. % Input :
  8. % A : the n times n matrix
  9. % b : the right hand side vector
  10. % tol : the stopping criterion
  11. % maxiter : the maximum number of iterations
  12. % x0 : the initial approximation
  13. % iprint : iprint = 1 to print some components of the
  14. % approximated solution and the norm of the residual
  15. % in the Matlab Command Window.
  16. % : iprint = 0 if we do not want to print the components
  17. % and the norm of the residuals in the matlab Command
  18. % Window.
  19. %
  20. % Output :
  21. % x : the solution
  22. % tnorm : the pseudo residual norm array
  23. % flag : the stopping criterion uses
  24. % if flag = 0 : we use tol
  25. % if flag = 1 : we use maxiter
  26. %
  27. n=length(A(:,1));
  28. if (nargin < 5)
  29. x = x0;
  30. else
  31. x = zeros(n,1);
  32. end
  33. if (nargin < 4)
  34. maxiter = n;
  35. end
  36. if (nargin < 3)
  37. tol = 1.d-7;
  38. end
  39. rs = zeros(n,1); s = zeros(n,1); c = zeros(n,1);
  40. p = 1:n; b = b-A*x;
  41. [~,indice] = max(abs(b(1:n)));
  42. beta = b(indice); D_inv = 1/beta; rs(1) = beta;
  43. dnorm = abs(beta); tnorm = dnorm;
  44. if ((indice < 1)||(indice > 1))
  45. temp = p(1); p(1) = p(indice); p(indice) = temp;
  46. temp = b(1); b(1) = b(indice); b(indice) = temp;
  47. Ai = A(:,indice);
  48. A(:,indice) = A(:,1); A(:,1) = Ai;
  49. Aip = A(indice,:); A(indice,:)=A(1,:); A(1,:) = Aip;
  50. end
  51. b = D_inv.*b;
  52. k = 0;
  53. while((dnorm > tol)&&(k < maxiter))
  54. k = k+1;
  55. n_k = n-k;
  56. k1 = k+1;
  57. if (iprint==1)
  58. fprintf(' %3.0f %6.4e \n',k-1,dnorm);
  59. end
  60. work = A(:,k);
  61. work = A(:,k1:n)*b(k1:n)+work;
  62. A(k1:n,k) = b(k1:n);
  63. for j = 1:k
  64. A(j,k) = work(j);
  65. work(j) = 0.d0;
  66. j1 = j+1;
  67. work(j1:n) = -A(j,k).*A(j1:n,j)+work(j1:n);
  68. end
  69. if (k < n)
  70. [~,indice] = max(abs(work(k1:n)));
  71. indice = k+indice;
  72. H_k1 = work(indice);
  73. if ((indice < k1)||(indice > k1))
  74. temp = p(k1); p(k1) = p(indice); p(indice) = temp;
  75. temp = work(k1); work(k1) = work(indice); work(indice) = temp;
  76. Ai = A(:,indice); A(:,indice) = A(:,k1);
  77. A(:,k1) = Ai; Aip = A(indice,:); A(indice,:) = A(k1,:);
  78. A(k1,:) = Aip;
  79. end
  80. if ((H_k1 < 0)||(H_k1 > 0))
  81. D_inv = 1/H_k1; b = work; b = D_inv.*b;
  82. end
  83. else
  84. H_k1 = 0;
  85. b = zeros(n,1);
  86. end
  87. % solving the least squares problem
  88. if (k > 1)
  89. for i = 2:k
  90. i_1 = i-1;
  91. dtemp = A(i_1,k);
  92. A(i_1,k) = c(i_1)*dtemp+s(i_1)*A(i,k);
  93. A(i,k) = -s(i_1)*dtemp+c(i_1)*A(i,k);
  94. end
  95. end
  96. if (k <= n)
  97. dgam = sqrt(A(k,k)^2+H_k1^2);
  98. c(k) = abs(A(k,k))/dgam;
  99. s(k) = (H_k1*abs(A(k,k)))/(A(k,k)*dgam);
  100. rs(k1) = -s(k)*rs(k);
  101. rs(k) = c(k)*rs(k);
  102. A(k,k) = c(k)*A(k,k)+s(k)*H_k1;
  103. dnorm = abs(rs(k1));
  104. end
  105. tnorm = [tnorm dnorm];
  106. end
  107. if (k==maxiter)
  108. flag=1;
  109. else
  110. flag=0;
  111. end
  112. if (iprint==1)
  113. fprintf(' %3.0f %6.4e \n',k,dnorm);
  114. end
  115. % computing d the solution of the least square problem
  116. for j = k:-1:1
  117. rs(j) = rs(j)/A(j,j);
  118. rs(1:j-1) = rs(1:j-1)-rs(j).*A(1:j-1,j);
  119. end
  120. work(1:k) = rs(1:k);
  121. for j = k:-1:1
  122. rs(j+1:k) = rs(j+1:k)+rs(j).*A(j+1:k,j);
  123. end
  124. c = A(k1:n,1:k)*work(1:k);
  125. x(1:k) = rs(1:k)+x(1:k);
  126. x(k1:n) = c(1:n_k)+x(k1:n);
  127. % Reorder the components of x to get the solution
  128. xx(p(1:n)) = x(1:n);
  129. x = xx';