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 - function x = lin_solve(A,b)
-
[R,v] = householder(A);
-
y = Qt_times_b(v,b);
-
x = R\y;
Here are the individual functions: - function [R,v] = householder(A)
-
[m,n] = size(A);
-
if m>=n,
-
NumberOfReflections = n;
-
else
-
NumberOfReflections = m - 1;
-
end
-
R = A;
-
v = cell(NumberOfReflections,1);
-
for k = 1:NumberOfReflections,
-
x = R(k:m,k);
-
xnorm = norm(x);
-
if xnorm>0,
-
% Compute the normal vector of the reflector
-
v{k} = -x;
-
v{k}(1) = v{k}(1) - sign(x(1))*xnorm;
-
v{k} = (sqrt(2)/norm(v{k}))*v{k};
-
% Update R
-
for j = k:NumberOfReflections,
-
R(k:m,j) = R(k:m,j) - (v{k}'*R(k:m,j))*v{k};
-
end
-
else
-
v{k} = zeros(m-k+1,1);
-
end
-
end
-
-
function y = Qt_times_b(v,b)
-
NumberOfReflections = length(v);
-
y = b;
-
p = NumberOfReflections;
-
for k = 1:NumberOfReflections,
-
F = eye(p)-2*(v{k}*v{k}')/norm(v{k})^2;
-
p = p -1;
-
% Put F into "Q" to get Qk
-
y = Qk*y;
-
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
0 1989 Sign in to post your reply or Sign up for a free account.
Similar topics
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...
|
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...
|
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
|
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...
|
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...
|
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...
|
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...
|
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...
|
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...
|
by: DolphinDB |
last post by:
The formulas of 101 quantitative trading alphas used by WorldQuant were presented in the paper 101 Formulaic Alphas. However, some formulas are complex, leading to challenges in calculation.
Take...
|
by: ryjfgjl |
last post by:
ExcelToDatabase: batch import excel into database automatically...
|
by: Vimpel783 |
last post by:
Hello!
Guys, I found this code on the Internet, but I need to modify it a little. It works well, the problem is this: Data is sent from only one cell, in this case B5, but it is necessary that data...
|
by: jfyes |
last post by:
As a hardware engineer, after seeing that CEIWEI recently released a new tool for Modbus RTU Over TCP/UDP filtering and monitoring, I actively went to its official website to take a look. It turned...
|
by: PapaRatzi |
last post by:
Hello,
I am teaching myself MS Access forms design and Visual Basic. I've created a table to capture a list of Top 30 singles and forms to capture new entries. The final step is a form (unbound)...
|
by: CloudSolutions |
last post by:
Introduction:
For many beginners and individual users, requiring a credit card and email registration may pose a barrier when starting to use cloud servers. However, some cloud server providers now...
|
by: Defcon1945 |
last post by:
I'm trying to learn Python using Pycharm but import shutil doesn't work
|
by: Shællîpôpï 09 |
last post by:
If u are using a keypad phone, how do u turn on JavaScript, to access features like WhatsApp, Facebook, Instagram....
|
by: af34tf |
last post by:
Hi Guys, I have a domain whose name is BytesLimited.com, and I want to sell it. Does anyone know about platforms that allow me to list my domain in auction for free. Thank you
| |