Arguments

function.If A is a function which returns A*x, it must
have the following header :

functiony=A(x)

If A is a function which returns A*x or A'*x depending t.
If t = "notransp", the function returns A*x.
If t = "transp", the function returns A'*x. It must
have the following header :

functiony=A(x, t)

Ap

function returning A'*x. It must have the followinf header :

functiony=Ap(x)

b

right hand side vector

x0

initial guess vector (default: zeros(n,1))

M1

left preconditioner : matrix or function (In the first case, default: eye(n,n)). If M1 is a function, she returns either,

only M1*x

or

M1*x or M1'*x depending t.

M1p

must only be provided when M1 is a function returning M1*x.
In this case M1p is the function which returns M1'*x.

M2

right preconditioner : matrix or function (In the first case, default: eye(n,n)). If M2 is a function, she returns either

only M2*x

or

M2*x or M2'*x depending t.

M2p

must only be provided when M2 is a function returning M2*x.
In this case M2p is the function which returns M2'*x

maxi

maximum number of iterations (default: n)

tol

error tolerance (default: 1000*%eps)

x

solution vector

flag

0 =

gmres converged to the desired tolerance within maxi iterations

1 =

no convergence given maxi

res

residual vector

err

final residual norm

iter

number of iterations performed

Description

Solves the linear system Ax=b using the Quasi Minimal Residual Method with preconditioning.

Examples

// If A is a matrixA=[94000028003200591350001000013723420000650534114000005500207002832120280000872003300000282071390001000320394680320001233088211006555000011100];b=ones(10,1);[x,flag,err,iter,res]=qmr(A,b)[x,flag,err,iter,res]=qmr(A,b,zeros(10,1),eye(10,10),eye(10,10),10,1d-12)// If A is a functionfunctiony=Atimesx(x, t)A=[94000028003200591350001000013723420000650534114000005500207002832120280000872003300000282071390001000320394680320001233088211006555000011100];if(t=='notransp')theny=A*x;elseif(t=='transp')theny=A'*x;endendfunction[x,flag,err,iter,res]=qmr(Atimesx,b)[x,flag,err,iter,res]=qmr(Atimesx,b,zeros(10,1),eye(10,10),eye(10,10),10,1d-12)// ORfunctiony=funA(x)A=[94000028003200591350001000013723420000650534114000005500207002832120280000872003300000282071390001000320394680320001233088211006555000011100];y=A*xendfunctionfunctiony=funAp(x)A=[94000028003200591350001000013723420000650534114000005500207002832120280000872003300000282071390001000320394680320001233088211006555000011100];y=A'*xendfunction[x,flag,err,iter,res]=qmr(funA,funAp,b)[x,flag,err,iter,res]=qmr(funA,funAp,b,zeros(10,1),eye(10,10),eye(10,10),10,1d-12)// If A is a matrix, M1 and M2 are functionsfunctiony=M1timesx(x, t)M1=eye(10,10);if(t=="notransp")theny=M1*x;elseif(t=="transp")theny=M1'*x;endendfunctionfunctiony=M2timesx(x, t)M2=eye(10,10);if(t=="notransp")theny=M2*x;elseif(t=="transp")theny=M2'*x;endendfunction[x,flag,err,iter,res]=qmr(A,b,zeros(10,1),M1timesx,M2timesx,10,1d-12)// ORfunctiony=funM1(x)M1=eye(10,10);y=M1*x;endfunctionfunctiony=funM1p(x)M1=eye(10,10);y=M1'*x;endfunctionfunctiony=funM2(x)M2=eye(10,10);y=M2*x;endfunctionfunctiony=funM2p(x)M2=eye(10,10);y=M2'*x;endfunction[x,flag,err,iter,res]=qmr(A,b,zeros(10,1),funM1,funM1p,funM2,funM2p,10,1d-12)// If A, M1, M2 are functions[x,flag,err,iter,res]=qmr(funA,funAp,b,zeros(10,1),funM1,funM1p,funM2,funM2p,10,1d-12)[x,flag,err,iter,res]=qmr(Atimesx,b,zeros(10,1),M1timesx,M2timesx,10,1d-12)