Sistemas Lineares com Múltiplos Lados Direitos e Operadores de
Projeção
Valéria S. Motta∗
Departamento Básico,
Instituto Militar de Engenharia - IME,
Praça General Tibúrcio, 80, Praia Vermelha,
22290-270, Rio de Janeiro, RJ.
E-mail: [email protected].
Luiz M. Carvalho
Departamento de Matemática Aplicada,
Instituto de Matemática e Estatística,
Universidade do Estado do Rio de Janeiro - UERJ,
Rua São Francisco de Xavier, 524, sala 6026, bloco D, Maracanã,
20550-013, Rio de Janeiro, RJ.
E-mail: [email protected].
Nelson Maculan
Programa de Engenharia de Sistemas e Computação - COPPE,
Universidade Federal do Rio de Janeiro - UFRJ,
C. P. 68511, 21945-970, Rio de Janeiro - RJ.
E-mail: [email protected].
1
Introdução
A resolução de sistemas lineares com grande número de incógnitas é uma área de muito interesse na Ciência da Computação.
Vários métodos foram desenvolvidos nessa
busca, podemos citar alguns deles: Lanczos
[?], [?], Gradiente Conjugado (CG) [?], [?],
Mínimo Resíduo Generalizado (GMRES) [?],
Gradiente Bi-Conjugado (Bi-CG) [?], Gradiente Conjugado Quadrado (CGS) [?], Gradiente Bi-Conjugado Estabilizado (Bi-CGSTAB)
[?], Gradiente Bi-Conjugado Estabilizado em
Bloco (Bi-CGSTAB em Bloco) [?].
A escolha do melhor método para se resolver um determinado sistema depende de algumas características da matriz principal do sistema, se ela é simétrica, positiva-definida, nãosimétrica, semi-definida, esparsa, entre outras.
Esses métodos iterativos são baseados em projeções em subespaços de Krylov e precisam
∗
Aluna de doutorado do Programa de Engenharia de Sistemas e Computação - COPPE - UFRJ.
Email:[email protected]
construir bases para esses espaços. Tais bases devem ser numericamente estáveis, então,
devem ser ortogonalizadas. A ortogonalização
de uma base de um subespaço de Krylov pode
ser conseguida através de, principalmente, dois
métodos básicos, o de Arnoldi [?], [?] e o da
Bi-Ortogonalização de Lanczos [?], [?].
Alguns dos algoritmos acima apresentam sua
versão em Bloco, a saber: o algoritmo de Lanczos em Bloco [?], o Gradiente Conjugado em
Bloco (BCG) (em mais de uma versão)[?], [?],
o Gradiente Bi-Conjugado em bloco (BBi-CG)
[?] .
Neste trabalho, as letras maiúsculas denotarão
matrizes; o índice superior em uma matriz, vetor ou em um escalar denota o número da iteração considerada; colunas de uma matriz são
indicadas pelo índice baixo; elementos de uma
matriz têm indicado no índice baixo o número
(k)
de sua linha e coluna. Assim, por exemplo xi,j
é o elemento da matriz X (k) da i-ésima linha
e j-ésima coluna. As letras α, β, γ, ε, η, ν, e ω
denotam matrizes de ordem s × s. A letra σi
denota os valores singulares, enquanto que λi
denota os autovalores. Em relação aos teoremas enunciados, no decorrer do trabalho, serão
indicadas referências para cada uma de suas demonstrações.
Esse trabalho é baseado principalmente em [?],
analisando os métodos a partir de operadores
de projeção, como contribuição algumas proposições novas são propostas e demonstradas.
Serão chamados de Teoremas os resultados que
se encontram em [?] e Proposições os resultados
que obtivemos.
Na segunda seção, faremos uma revisão sobre
operadores de projeção, nessa seção apresentamos a primeira proposição.
Na terceira seção definimos um subespaço de
Krylov e o método e algoritmo de Arnoldi [?],
[?], a segunda proposição é apresentada e sob
seu enfoque o algoritmo é reescrito.
Na quarta seção abordamos os algoritmos de
Lanczos [?], [?], teoremas e conceitos relacionados, segue daí a terceira proposição e o algoritmo de Lanczos Bi-Ortogonal [?], [?], [?]
reescrito.
Na quinta seção, o algoritmo do Gradiente BiConjugado Estabilizado em Bloco [?], [?], [?] é
o exemplo de algoritmo para resolução de sistemas com múltiplos lados direitos. Esse algoritmo é estudado e a partir dele temos a quarta
proposição.
2
Operadores de Projeção
Px ∈ K
x
x − P x⊥L
L
K
Px
Figura 1: Projeção de x sobre M e ortogonal a
L
bespaços resulta em Cn , então, existe uma única
projeção P onde Im(P ) = M e N (P ) = S.
Prova [?].
Essa transformação linear P é dita projeção sobre M , ao longo ou paralela ao
subespaço linear S.
O espaço vetorial
Cn se expressa como a soma direta da
Im(P ) = M e N (P ) = S, sabemos que
N (P ) = Im(I − P ). Se dim(Im(P )) = m
então, dim(N (P )) = dim(Im(I − P )) = n − m.
Portanto, podemos definir S a partir de
seu complemento ortogonal L = S ⊥ , onde
dim(L) = dim(S ⊥ ) = m.
(
P (x) ∈ M
(3)
x − P (x)⊥L.
As equações dadas em (??) definem uma projeção sobre M e ortogonal ao subespaço L. Uma
Definição 1 Considere P uma transformação projeção genérica é ilustrada na Figura (1).
linear, P é dita uma projeção, P : Cn → Cn se
P 2 = P , P é chamado também idempotente.
2.0.1 Métodos Gerais de Projeção
Todo elemento de Cn (espaço vetorial complexo
Sejam A uma matriz real de ordem n × n e,
n-dimensional) pode ser escrito na forma,
ainda, K e L dois subespaços m-dimensionais
n
x = P (x) + (I − P )(x)
(1) de R . Uma técnica de projeção sobre o subespaço K e ortogonal a L é um processo de
portanto, pode-se provar que, o espaço Cn pode encontrar a solução aproximada x̃ de Ax = b
ser decomposto na seguinte soma direta:
tendo condições para que x̃ pertença a K e a de
que o novo vetor resíduo seja ortogonal a L.
Cn = N (P ) ⊕ Im(P ),
(2)
Determine x̃ ∈ K tal que b − Ax̃⊥L.
onde N (P ) = {x ∈ Cn |P (x) = 0} é o núcleo de
P e Im(P ) = {y ∈ Cn |∃x ∈ Cn ; P (x) = y}, a Se desejarmos utilizar uma solução inicial x0 ,
imagem de P .
então a solução aproximada deve ser procurada
A recíproca é dada no teorema a seguir.
no espaço afim x0 + K ao invés de no espaço
vetorial K. O problema pode ser redefinido por:
Teorema 1 Para todo par de subespaços M e
S de Cn , tais que, a soma direta desses dois su- Determine x̃ ∈ x0 + K tal que b − Ax̃⊥L.
Observe que, se x̃ é escrito na forma x̃ = x0 + δ e
e o vetor resíduo inicial r0 é definido como
r0 = b − Ax0 ,
δ = V (W H AV )−1 W H Ae.
(10)
(4)
Considere P = V (W H AV )−1 W H A, podemos
2
então a equação acima é dada por b − A(x0 + verificar que P = P , ou seja, P é a projeção
dada na Proposição 1.
δ)⊥L ou
Da equação (??) segue que δ = P e, por hipór0 − Aδ⊥L.
(5) tese temos que δ ∈ K, então podemos concluir
que δ é a projeção do erro e sobre o subespaço
Em outras palavras, a solução aproximada pode K.⋄
Para que o Algoritmo 2.1 funcione, devemos ter
ser definida como:
que a matriz W H AV é invertível.
x̃ = x0 + δ, δ ∈ K
(6)
(r0 − Aδ, w) = 0, ∀w ∈ L.
(7) Teorema 2 Sejam A, K e L satisfazendo às
condições abaixo:
Esse é o passo básico da projeção, na sua
forma mais geral. A maioria das técnicas usam
uma sucessão de projeções. Normalmente, cada
novo passo da projeção usa um novo par de subespaços K e L e uma nova solução inicial x0 ,
essa, igual a mais recente aproximação obtida
no passo anterior da projeção. Métodos projetivos formam uma estrutura para os principais métodos iterativos utilizados na resolução
se sistemas lineares.
Um protótipo da técnica de projeção é representada no algoritmo a seguir.
Algoritmo 2.1. Protótipo do Método de
Projeção
1. Sejam a matriz A e os vetores x0 e b, Do
2.
Selecione um par de subespaços K e L
3.
Escolha bases V = [v1 , . . . , vm ] e
W = [w1 , . . . , wm ] para K e L
4.
r := b − Ax
5.
y := (W H AV )−1 W H r
6.
x := x + V y
7. EndDo
(i) A é positiva definida e K = L ou
(ii) A é não-singular e L = AK.
Então a matriz W H AV é não-singular para
quaisquer bases de V e W de K e L, respectivamente.
Prova: [?].
3
3.1
Método de Arnoldi
Subespaço de Krylov
Definição 2 O espaço m-dimensional gerado
pelo vetor v e pelas potências de A multiplicadas por v até o expoente (m−1), é chamado subespaço de Krylov m-dimensional [?], [?], [?],
notamos tal subespaço por Km (A, v).
Os métodos iterativos que serão abordados buscam a solução aproximada do sistema Ax = b
em um subespaço de Krylov conveniente. Considerando x0 = 0 a solução inicial, xi a solução
aproximada e ri = b − Axi o resíduo, ambos na
Proposição 1 Seja o Algoritmo 2.1. Então, iteração i, sabemos que:
existe um operador de projeção P , que projeta
ri = pi (A)r0 .
(11)
o erro e sobre o subespaço K. A matriz que
representa tal operador é dada por:
onde pi (A) é um polinômio de grau i. Segue
H
−1
H
P = V (W AV ) W A
(8) que:
i
X
Prova:
(I − A)j r0 ∈ [r0 , Ar0 , . . . , Ai r0 ] ≡
xi+1 =
O passo (6) do algoritmo é semelhante à equaj=0
ção (??), tomemos δ = V y.
≡ Ki+1 (A, r0 ).
(12)
No passo (7) temos y = (W H AV )−1 W H r, sabendo que r = Ae, onde e é o erro, temos que: Nessa abordagem o subespaço de Krylov é o
conjunto ao qual pertencem as soluções aproxiH
−1
H
y = (W AV ) W Ae,
(9) madas de um sistema linear da forma, Ax = b.
Vamos obter uma base para esse subespaço.
Prova:
Uma base que deve ser obviamente considerada Consideremos V uma base para o subespaço de
é:
Krylov, Km (A, v1 ). Considerando do passo (2)
ao passo (7), teremos:
{r0 , Ar0 , . . . , Ai−1 r0 } ∈ Ki (A, r0 )
(13)
j
X
hi,j vi
(16)
w
=
Av
−
mas essa base pode não ser muito atrativa sob
j
i=1
o ponto de vista numérico. Os vetores A r0
Pj
apontam mais e mais para a direção do autoFaçamos Av = x e P (x) =
i=1 hi,j vi , mas
vetor associado ao maior autovalor, quando j
P (x) pode ser escrita como P (x) = V y. Quecresce. Com isso, a base dos vetores pode ir
remos provar que P é uma projeção, ou seja,
se tornando linearmente dependente em aritP 2 = P . Temos que:
mética finita.
Na próxima seção veremos como construir uma
w = x−Vy
(17)
base adequada para um subespaço de Krylov.
H
H
H
V w = V x−V Vy
Mas , por hipótese, as colunas de V são ortogonais a w, assim V H w = 0 e temos também que
O processo de Arnoldi é um método para se V H V = I, então:
computar uma base ortonormal {v1 , v2 , . . . , vj }
do subespaço de Krylov, Kj (A, v1 ). A ortonory = V H x, como P (x) = V y,
malização admite diferentes procedimentos, um
temos P = V V H ,
deles é o Gram-Schmidt modificado.
observe que P 2 = P , logo P é a projeção dada
na Proposição 2.
Algoritmo 2.2. Algoritmo de Arnoldi
Podemos concluir que a projeção P = V V H
1. v1 = krr00k2
projeta o vetor Avi sobre o subespaço de
2. f or j = 1, . . . , m − 1
Krylov, Ki (A, v1 ), com i = 1, . . . , m.⋄
3.
wj+1 = Avj ;
3.2
Método de Arnoldi
4.
f or i = 1, . . . , j
5.
hi,j = viT wj+1 ;
6.
wj+1 = wj+1 − hi,j vi ;
7.
end;
8.
hj+1,j = kwj+1 k2 ;
wj+1
;
9.
vj+1 = hj+1,j
10. end.
4
Métodos de Lanczos
O algoritmo simétrico de Lanczos [?], [?], [?],
pode ser visto como uma simplificação do método de Arnoldi para o caso particular em que
Em termos matriciais, os passos (2) a (7) do a matriz considerada é simétrica. Quando A é
simétrica, então a matriz de Hessenberg Hm é
Algoritmo 2.2 podem ser escritos como:
uma tridiagonal simétrica.
VmH AVm = Hm
(14)
onde Hm é uma matriz de Hessenberg superior,
ou seja, são nulos todos os elementos abaixo da
primeira subdiagonal (diagonal logo abaixo da
diagonal principal).
Algoritmo 2.3. Algoritmo de Lanczos
1. Escolha um vetor v1 tal que kv1 k2 = 1. Seja
β1 ≡ 0, v0 ≡ 0.
2. Para j = 1, 2, . . . , m, Do
3.
wj = Avj − βj vj−1
4.
αj = (wj , vj )
5.
wj = wj − αj vj
6.
βj+1 = kwj k2 . Se βj+1 = 0 então Stop
wj
7.
vj+1 = βj+1
8. EndDo
Proposição 2 Seja o Algoritmo 2.2. Então,
existe um operador de projeção P , que projeta o vetor Avi sobre o subespaço de Krylov
Ki (A, v1 ), onde i = 1, . . . , m . A matriz que
representa tal operador é dada por:
Relacionando o método de Lanczos com a projeção, os resultados são análogos ao método de
P =VVH
(15) Arnoldi, sua projeção é dada por P = V V H e
projeta o vetor Avi sobre o subespaço de Krylov dado por Ki (A, v1 ), com i = 1, . . . , m. A
particularidade é que a matriz H é tridiagonal,
como já vimos anteriormente.
Em aritmética exata, o ponto central do algoritmo de Lanczos é a relação da forma:
Nessa formulação cada vetor das duas bases, V
e W , têm norma unitária e (vj , wj ) = 1, mas
os vetores podem ser também escolhidos de tal
forma que (vj , wj ) = dj .
Teorema 3 Se não ocorre falha no Algoritmo
(18) (2.4) antes de m passos, então:
(
Nessa recorrência temos três vetores envolvidos,
1, se i = j,
vj−1 , vj e vj+1 , por isso é chamada recorrência (vi , wj ) =
, 1 ≤ i, j ≤ m. (19)
0, se i 6= j.
de três termos [?], [?], [?]. Temos que guarβj+1 vj+1 = Avj − αj vj − βj vj−1 .
dar três vetores a cada iteração no Algoritmo
de Lanczos, enquanto que no Algoritmo de Ar- Prova: [?].
noldi guardamos todos os vetores da base do A condição de bi-ortogonalidade implica que:
subespaço de Krylov calculados até a iteração
V HW = I
(20)
corrente.
Vamos considerar agora a matriz A nãoProposição 3 Seja o Algoritmo 2.4. Então,
simétrica.
existe um operador de projeção P , que projeta
o vetor Avi sobre o subespaço de Krylov
Definição 3 Dois conjuntos de vetores U =
{u1 , . . . , um } e V = {v1 , . . . , vm } são ditos bi- Ki (A, v1 ), onde i = 1, . . . , m . A matriz que
representa tal operador é dada por:
ortogonais se:
(ui , vj ) = δij , onde
(
1, se i = j,
deltaij =
0, se i 6= j.
Na forma matricial temos V H U = I.
Não poderemos resolver o sistema Ax = b
com uma única recorrência de três termos,
mas podemos obter uma base adequada, nãoortogonal, com uma recorrência de três termos,
desde que essa base seja bi-ortogonal em relação a alguma outra base dada anteriormente.
Vamos construir os conjuntos de vetores,
{v1 , . . . , vi } ⊂ Ki (A, r0 ) e {w1 , . . . , wi } ⊂
Ki (AH , r̂0 ), partindo de um dos algoritmos de
Lanczos(pag. 34 [?]).
Algoritmo 2.4. Algoritmo da
Bi-Ortogonalização de Lanczos
1. Escolha dois vetores v1 , w1
tais que (v1 , w1 ) = 1
2. Sejam β1 ≡ δ1 ≡ 0, w0 ≡ v0 ≡ 0
3. Para j = 1, 2, . . . , m. Do
4.
αj = (Avj , wj )
5.
ṽj+1 = Avj − αj vj − βj vj−1
6.
w̃j+1 = AH wj − αj wj − δj wj−1
7.
δj+1 = |(ṽj+1 , w̃j+1 )|1/2 . If δj+1 = 0 Stop
8.
βi+1 = (ṽj+1 , w̃j+1 )/δj+1
9.
wj+1 = w̃j+1 /βj+1
10.
vj+1 = ṽj+1 /δj+1
11. EndDo
P = V WH
(21)
Temos também, o operador de projeção Q =
P H , que projeta o vetor AH wi sobre o subespaço
de Krylov Ki (AH , w1 ), onde i = 1, . . . , m . A
matriz que representa tal operador é dada por:
Q = WV H
(22)
Prova:
Partindo
do passo (5) e fazendo x = Avj e
P
α
v
−
βj vj−1 = P (x) = V y, temos:
j j j
x̃ = x − V y
(23)
W H x̃ = W H x − W H V y
(24)
façamos:
como as colunas de W são ortogonais a x̃ e
W H V = I, segue que:
y = W H x assim P (x) = V W H x.
(25)
Observe que P 2 = P , assim P é uma projeção e
sua matriz associada V W H projeta o vetor Avj
sobre o subespaço de Krylov Kj = K(A, v1 ),
com, j = 1, . . . , m.
Analogamente
podemos
mostrar
que
Q = W V H .⋄
5
Gradiente Bi-Conjugado Estabilizado em Bloco
Prova:
Observando os passos (6) e (7), podemos
reescrevê-los como:
Um dos objetivos dos algoritmos em bloco [?],
X̂ = X + ∆
(28)
[?], [?] é o de resolver sistemas com múltiplos
R̂ = R − A∆,
(29)
lados direitos de forma mais eficiente do que os
algoritmos que não são em bloco. Nem sempre isso é possível, mas quando a matriz A é com ∆ ∈ K(A, R0 ), tomamos ∆ = V Y .
densa, ou quando trabalhamos com um pré- Por hipótese o resíduo é ortogonal ao subespaço
condicionador o método pode ser bem atrativo de Krylov, L = (AH , R̃0 ), assim:
[?].
W H R̂ = W H R − W H AV Y
(30)
Vamos considerar o seguinte sistema:
AX = B,
(26) Mas W H R̂ = 0, assim:
Y = (W H AV )−1 W H R
(31)
onde A é uma matriz real de ordem n × n , nãosingular e não-simétrica, B e X são matrizes de
sabemos que o resíduo R é dado por R = AE,
ordem n × s.
O método iterativo toma uma solução X0 e o onde E é o erro, temos que:
resíduo inicial R0 = B − AX0 , com X0 , R0 ∈
Y = (W H AV )−1 W H AE,
(32)
Rn×s . No passo k teremos Xk = X0 + ∆k uma
aproximação da solução do sistema (??), onde
∆k ∈ Kk (A, R0 ) = [R0 , AR0 , . . . , Ak−1 R0 ] que e
é um subespaço de Krylov em Bloco. Vejamos
∆ = V (W H AV )−1 W H AE.
(33)
o algoritmo do Gradiente Bi-Conjugado Estabilizado em Bloco [?]:
Considere Q = V (W H AV )−1 W H A, podemos
verificar que Q2 = Q, ou seja, Q é uma projeAlgoritmo 2.5. Algoritmo do Gradiente
ção que projeta o erro ∆ sobre o subespaço de
Bi-Conjugado Estabilizado em Bloco
Krylov K = K(A, R0 ).⋄
1. Sejam X0 é uma solução inicial e
R0 = B − AX0 o resíduo, ambas matrizes n × s e
R0 de posto s.
2. Seja R̃0 uma matrix n × s arbitrária de posto
s.
3. P0 = R0 , P̃0 = R̃0 .
4. Para j = 0, 1, . . . Do
5.
αj = (P̃jH APj )−1 R̃jH Rj
6.
Xj+1 = Xj + Pj αj
7.
Rj+1 = Rj − APj αj
8.
α̃j = (PjH AH P̃j )−1 RjH R̃j
9.
R̃j+1 = R̃j − AH P̃j α̃j
H
10.
βj = (R̃jH Rj )−1 R̃j+1
Rj+1
H
−1 H
11.
β̃j = (Rj R̃j ) Rj+1 R̃j+1
12.
Pj+1 = Rj+1 + Pj βj
13.
P̃j+1 = R̃j+1 + P̃j β̃j
14. EndDo
Proposição 4 Seja o Algoritmo 2.5. Então,
existe um operador de projeção Q, que projeta
o erro E sobre o subespaço de Krylov, K =
K(A, R0 ). A matriz que representa tal operador
é dada por:
Q = V (W H AV )−1 W H A
(27)
6
Conclusões
Nesse trabalho apresentamos uma diferenciação
entre os métodos iterativos baseados no método
de Arnoldi e os que são baseados no método de
Lanczos Bi-Ortogonal.
Para tal, estudamos os algoritmos clássicos sob
o enfoque de projeções. Assim, alguns algoritmos foram apresentados de uma nova forma,
usando os operadores de projeção implícitos em
cada método abordado.
Temos como objetivos futuros implementar alguns métodos em bloco, na resolução de sistemas com múltiplos lados direitos e compará-los.
Referências
[1] H. Dai. Block Bidiagonalization Methods
For Solving Nonsymmetric Linear Systems
With Multiple Right-Hand Sides. Technical report, 1998.
[2] R. Fletcher. Conjugate Gradient Methods [13] J. Rice, editor. The block Lanczos method
for Indefinite Systems. In G. Watson, edifor computing elgenvalues, New York,
tor, Proc. of the Dundee Biennial Confe1977. Academic Press.
rence on Numerical Analysis (1975), volume 506 of Lecture Notes in Mathematics, [14] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2nd ed. edition, 2003.
pages 73–89. Springer Verlag„ Berlin, Heidelberg, New York, 1975.
[15] Y. Saad and M. H. Schultz. GMRES:
A generalized minimal residual algorithm
[3] G. H. Golub and D. P. O’Leary. Some Hisfor solving nonsymmetric linear systems.
tory of the Conjugate Gradient and LancSIAM Journal on Scientific and Statistical
zos Algorithms: 1948–1976. SIAM Review,
Computing, 7(3):856–869, 1986.
31(1):50–102, 1989.
[16] J. R. Shewchuk. An Introduction to the
Conjugate Gradient Method Without the
Agonizing Pain. Technical report, Pittsburgh, PA, USA, 1994.
[5] A. E. Guennouni and K. Jbilou. A Block
Version of BICGSTAB for Linear Systems [17] P. Sonneveld. CGS, A Fast Lanczos-Type
with Multiple Rigth -Hand Sides. ETNA,
Solver for Nonsymmetric Linear Systems.
16:129–142, 2002.
SIAM Journal on Scientific and Statistical
[4] A. Greenbaum. Iterative Methods for Solving Linear Systems. SIAM, 1997.
Computing, 10(1):36–52, 1989.
[6] A. E. Guennouni, K. Jbilou, and H. Sadok.
The block Lanczos method for [18] H. A. van der Vorst. BI-CGSTAB: A Fast
linear systems with multiple right-hand
and Smoothly Converging Variant of BIsides. Applied Numerical Mathematics,
CG for the solution of Nonsymmetric Li51(2-3):243–256, 2004.
near Systems. SIAM Journal on Scientific
and Statistical Computing, 13(2):631–644,
[7] M. R. Hestenes and E. Stiefel. Methods
1992.
of Conjugate Gradients for Solving Linear
Systems. Journal of Research os the Nati- [19] H. A. van der Vorst. Iterative Krylov
onal Bureau of Standards, 49(6):409–436,
Methods for a Large Systems. Cambridge
December 1952.
University Press, 2003.
[8] I. C. F. Ipsen and C. D. Meyer. The Idea
Behind Krylov Methods. American Mathematical Monthly, 105(10):889–899, 1998.
[9] C. Lanczos. An iterarion method for the
solution of the eigenvalue problem of linear
differential and integral operators. J. Res.
Nat. Bur. Standards, (45):255–282, 1950.
[10] Lanczos.C. Solution of systems of linear
equations by minimized iterations. J. Res.
Natl. Bur. Stand, 49:pp. 33–53, 1952.
[11] D. O’Leary. The block conjugate gradient algorithm and related methods. Linear
Algebra and Its Applications, 29:293–322,
1980.
[12] C. C. Paige and M. A. Saunders. Solution
of sparse indefinite systems of linear equations. SIAM J. Numer. Anal., 12:617–629,
1975.