473,425 Members | 1,856 Online
Bytes | Software Development & Data Engineering Community
Post Job

Home Posts Topics Members FAQ

Join Bytes to post your question to a community of 473,425 software developers and data experts.

Matlab: solving the linear system with QR/Householder

Hi all,
I'm trying to implement the QR method for solving the linear system Ax = b. The QR factorization is achieved using Householder method.

The main function is

Expand|Select|Wrap|Line Numbers
  1. function x = lin_solve(A,b)
  2. [R,v] = householder(A);
  3. y = Qt_times_b(v,b);
  4. x = R\y;
Here are the individual functions:

Expand|Select|Wrap|Line Numbers
  1. function [R,v] = householder(A)
  2. [m,n] = size(A);
  3. if m>=n,
  4.     NumberOfReflections = n;
  5. else
  6.     NumberOfReflections = m - 1;
  7. end
  8. R = A;
  9. v = cell(NumberOfReflections,1);
  10. for k = 1:NumberOfReflections,
  11.     x = R(k:m,k);
  12.     xnorm = norm(x);
  13.     if xnorm>0,
  14.         % Compute the normal vector of the reflector
  15.         v{k} = -x;
  16.         v{k}(1) = v{k}(1) - sign(x(1))*xnorm;
  17.         v{k} = (sqrt(2)/norm(v{k}))*v{k};
  18.         % Update R
  19.         for j = k:NumberOfReflections,
  20.                 R(k:m,j) = R(k:m,j) - (v{k}'*R(k:m,j))*v{k};
  21.         end
  22.     else
  23.         v{k} = zeros(m-k+1,1);
  24.     end
  25. end
  26.  
  27. function y = Qt_times_b(v,b)
  28. NumberOfReflections = length(v);
  29. y = b;
  30. p = NumberOfReflections;
  31. for k = 1:NumberOfReflections,
  32.    F = eye(p)-2*(v{k}*v{k}')/norm(v{k})^2;
  33.    p = p -1;
  34.    % Put F into "Q" to get Qk
  35.    y = Qk*y;
  36. end
I'm having trouble with the "Put F into Q to get Q_k" part.

I know Q is implicitly defined by the Householder reflectors stored in v, computed in the householder function, and y can be obtained by Q'b = Q_n * ... * Q_3 * Q_2 * Q_1 * b

Also, Qk =
[ I_(k-1) 0 ]
[ 0 F_k ]
matrix, where F_k = I - 2 (v * v') / norm(v)^2

But how do I put F into Q?

Thank you.

Regards,
Rayne
Sep 16 '07 #1
0 1994

Sign in to post your reply or Sign up for a free account.

Similar topics

2
by: Anatoly | last post by:
Hi All !!! Can anyone give me link to the C++ source or to the compiled program unit (compatible with Microsoft VC++), designed for the solving of the algebraic linear systems with rare matrices...
6
by: Lars Christiansen | last post by:
Hi! I am a master student in (geo)physics at the University of Copenhagen and part of a study group on C++ as a scientific programming language. I, and the other students in the group, have...
4
by: Thomas | last post by:
Hi, I've to solve a system of linear equations and inequations in C++. Are there (free) librarys for that topic? Regards, Thomas
0
by: m31hu1 | last post by:
Hello, i try to execute a .m file using PHP i wanted to count the number of jpg files in the directory but PHP was unable to detect it. the PHP script: <?php execute('start /b matlab...
6
by: omesh | last post by:
hi guys.. i got a problem.. could anyone help me write a program.. the program is as follows:- 2. Solving systems of Linear equations using Iteration: You are required to write a C++ program...
4
by: itcecsa | last post by:
Hi, I am implementing a small Python project, what I am going to do is to open Matlab and run some M-files, and get some output from Matlab command prompt. I have no idea how to open Matlab...
0
by: ursskmdali | last post by:
Hi, I have just created a sample in vc++ 6.0 to call MATLAB functions. Everything goes fine it compiles, builds fine with no errors and while executing it giving error like this "Missing ICU DATA...
5
by: adinda | last post by:
So what i need is this; (I'm very new at this,, programming in C I mean...) In matlab I had a while loop, and after each loop was done I added my resulting matrix to an object. Seeing the loop...
14
by: Luna Moon | last post by:
Dear all, Can C++/STL/Boost do the vectorized calculation as those in Matlab? For example, in the following code, what I really want to do is to send in a vector of u's. All other...
1
by: nemocccc | last post by:
hello, everyone, I want to develop a software for my android phone for daily needs, any suggestions?
0
marktang
by: marktang | last post by:
ONU (Optical Network Unit) is one of the key components for providing high-speed Internet services. Its primary function is to act as an endpoint device located at the user's premises. However,...
0
by: Hystou | last post by:
Most computers default to English, but sometimes we require a different language, especially when relocating. Forgot to request a specific language before your computer shipped? No problem! You can...
0
Oralloy
by: Oralloy | last post by:
Hello folks, I am unable to find appropriate documentation on the type promotion of bit-fields when using the generalised comparison operator "<=>". The problem is that using the GNU compilers,...
0
tracyyun
by: tracyyun | last post by:
Dear forum friends, With the development of smart home technology, a variety of wireless communication protocols have appeared on the market, such as Zigbee, Z-Wave, Wi-Fi, Bluetooth, etc. Each...
0
agi2029
by: agi2029 | last post by:
Let's talk about the concept of autonomous AI software engineers and no-code agents. These AIs are designed to manage the entire lifecycle of a software development project—planning, coding, testing,...
0
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 1 May 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome a new...
0
by: conductexam | last post by:
I have .net C# application in which I am extracting data from word file and save it in database particularly. To store word all data as it is I am converting the whole word file firstly in HTML and...
0
by: adsilva | last post by:
A Windows Forms form does not have the event Unload, like VB6. What one acts like?

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.