%% Antes de processar este arquivo LaTeX (LaTeX2e) deve
%% verificar que o arquivo TEMA.cls estah no mesmo
%% diretorio. O arquivo TEMA.cls pode ser obtido do
%% endereco www.sbmac.org.br/tema.

\documentclass{TEMA}

\usepackage[brazil]{babel}      % para texto em Português
%\usepackage[english]{babel}    % para texto em Inglês

\usepackage[latin1]{inputenc}   % para acentuação em Português
%\input{P-margin.inf}

\usepackage[dvips]{graphics}
\usepackage{subfigure}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{multicol}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage[centertags]{amsmath}
\usepackage{graphicx,indentfirst,amsmath,amsfonts,amssymb,amsthm,newlfont}
\usepackage{longtable}
\usepackage{cite}
\usepackage[usenames,dvipsnames]{color}
\usepackage{scalefnt}
\usepackage{float}



\newcommand{\B}{{\tt\symbol{92}}}
\newcommand{\til}{{\tt\symbol{126}}}
\newcommand{\chap}{{\tt\symbol{94}}}
\newcommand{\agud}{{\tt\symbol{13}}}
\newcommand{\crav}{{\tt\symbol{18}}}

\begin{document}

%********************************************************
\title
    {Sistemas lineares aproximados derivados de problemas de fluxo multiproduto em métodos de Pontos Interiores\thanks{Os autores agradecem à FAPESP e ao CNPq pelo suporte financeiro deste trabalho.}
}

\criartitulo

\runningheads{Autores}{Sistemas lineares aproximados em métodos de Pontos Interiores}

\begin{abstract}
{\bf Resumo}. Uma das abordagens utilizadas para resolver o sistema linear que surge a cada iteração nos métodos de pontos interiores primal-dual é reduzi-lo a um sistema linear equivalente simétrico definido positivo, conhecido como sistema de equações normais, e aplicar a fatoração de Cholesky na matriz do sistema. 

A grande  desvantangem desta abordagem é o preenchimento gerado durante a fatoração, o que pode tornar seu uso inviável, por limitação de tempo e memória. Com o intuito de contornar o problema de preenchimento gerado na fatoração de Cholesky, neste trabalho, estamos propondo uma abordagem que resolve de forma direta sistemas lineares aproximados do sistema de equações normais derivados de problemas de fluxo multiproduto e que exerce um certo controle sobre o preenchimento. 


{\bf Palavras-chave}.  Método de pontos interiores primal-dual,
Fatoração de Cholesky, Sistema de Equações Normais.
\end{abstract}


%********************************************************
\newsec{Introdução}

Um Problema de Programação Linear (PPL) consiste em encontrar uma solução que satisfaz  equações e/ou inequações lineares e que simultaneamente maximize ou minimize uma função linear denominada função objetivo.

Em sua forma padrão este problema é dado por
\begin{equation}
\begin{array}{lc} 
\text{min} & c^Tx\\
 \text{s.a} & Ax=b \\
  & x \ge 0, \\ 
 \end{array}\label{eqprimal}
\end{equation}
\noindent em que $A$ é uma matriz de ordem $m\times n$, $x\in\mathbb{R}^{n}$, $b\in\mathbb{R}^{m}$, $c\in\mathbb{R}^n$ e $c^Tx$ é a função objetivo do problema que deve ser minimizada.

A todo PPL na forma padrão, pode-se associar um problema dual que possui a seguinte forma

\begin{equation}
\begin{array}{lc} 
\text{max} &  b^{T}y\\
 & A^{T}y+z= c \\
  & y\in\mathbb{R}^{m}, z \geq 0, \\
 \end{array}\label{eqdual}
\end{equation} 
em que $z\in\mathbb{R}^{n}$ e $b^Ty$ é a função objetivo do problema dual que deve ser maximizada.
 
Por estarem associados desta forma, denominamos o par de PPLs (\ref{eqprimal}) e (\ref{eqdual}) de par primal-dual, sendo (\ref{eqprimal}) o problema primal e (\ref{eqdual}) o problema dual. 

Uma solução factível para o par primal-dual é um vetor $(x,y,z)$ que satisfaz as retrições  $Ax=b,$ $x\geq 0\,$ $A^{T}y+z= c$ e $z\geq 0$.  Uma solução ótima $(x^*, y^*, z^*)$ para o par primal-dual é uma solução factível com o menor e maior valor possível para as funções objetivos dos problemas primal e dual respectivamente. 

O teorema \ref{teo_otim}, enunciado a seguir, apresenta as condições de Karush-Kuhn-Tucker (KKT) para problemas de progração linear, que oferece uma condição necessária e suficiente para que um vetor $(x,y,z)$ seja uma solução ótima para o par primal-dual.

\begin{teoTEMA}\label{teo_otim}\textbf{(Condições de otimalidade)}\cite{WR96} $(x,y,z)$  são soluções ótimas para o par primal-dual se, e somente se, as seguintes condições são satisfeitas
\begin{align}
Ax=b,& \label{KKT1}\\
  A^{T}y+z= c, & \label{KKT2}\\
   x_iz_i=0,  & \hspace{0.3cm}  i=1,...,n, \label{KKT3}\\
   (x,z) \geq 0. \label{KKT4}
\end{align}
\end{teoTEMA}

\subsection{Método de Pontos Interiores Primal-Dual}

Os métodos de pontos interiores (MPI) para problemas de programação linear, como a maioria dos algoritmos iterativos em otimização, possuem dois pontos-chave: um procedimento que determina uma direção do passo ou direção de busca em cada iteração e uma medida conveniente deste passo. 

Nos métodos de pontos interiores do tipo primal-dual a direção de busca é calculada por meio de um sistema linear obtido  pela aplicação de variações do método de Newton \cite{Va96} às três igualdades (\ref{KKT1}), (\ref{KKT2}) e (\ref{KKT3}) \cite{WR96}.

Logo, a direção de busca $(dx, dy, dz)$ em cada iteração pode ser calculada por meio do sistema linear
\begin{equation}\label{eq1}
\begin{bmatrix}
A & 0 & 0 \\ 
0 & A^{t} & I  \\ 
Z & 0 & X%
\end{bmatrix}%
\begin{bmatrix}
dx \\ 
dy \\ 
dz%
\end{bmatrix}%
=%
\begin{bmatrix}
r_{p} \\ 
r_{d} \\ 
r_{c}%
\end{bmatrix}% 
\end{equation}
\noindent em que $X$ e $Z$ são matrizes diagonais, cujos elementos da diagonal são dados ordenadamente pelas entradas dos vetores $x$ e $z$, repsectivamente,  os resíduos são da forma $r_p=b-Ax$, $r_{d}=c-A^t y-z$ e $r_{c}=-XZe+\sigma\mu e$, com $\mu=x^t z/n$ (\textit{gap} de dualidade) e $\sigma$ é um parâmetro de centralização entre $0$ e $1$. 

O tamanho do passo é dado de forma que as variáveis $x$ e $z$ sejam estritamente positivas.

Em aplicações reais o sistema linear (\ref{eq1}) quase sempre possui dimensões elevadas e um alto grau de esparsidade e sua resolução é a etapa do método que requer mais esforço computacional. 

Por meio de operações algébricas, o sistema  pode ser reduzido a um sistema linear equivalente cuja matriz é simétrica definida positiva dado por 

\begin{equation}
ADA^Tdy=h. \label{eqnormais}
\end{equation}

em que $D=XZ^{-1}$ e
\begin{equation}
h=r_p+ADr_d-ADX^{-1}r_c.\label{h}
\end{equation} 

Esse sistema linear é conhecido como sistema de equações normais.
Uma das abordagens utilizadas para resolvê-lo é a aplicação da fatoração de Cholesky na matriz dos coeficientes, já que ela é simétrica definida positiva. 

A fatoração de Cholesky de uma matriz simétrica definida positiva $\mathcal{A}$ de ordem $m$ consiste em calcular uma matriz $L$ triangular inferior, denominada fator de Cholesky, tal que $\mathcal{A} =LL^T$. Os elementos da diagonal de $L$ são denominados pivôs.

 A  solução de um sistema linear com a matriz $\mathcal{A}$, pode então, ser calculada por  meio da resolução de dois sistemas triangulares com as matrizes $L$ e $L^T$. Esta abordagem é uma forma bastante estável para calcular soluções de sistemas lineares com a matriz dos coeficientes simétrica definida positiva, por isso, ela é utilizada nos MPI primal-dual sempre que a decomposição não é extremamente cara. 

A grande desvantagem da fatoração de Choleky é o preenchimento gerado durante a decomposição, que consiste no surgimento de elementos não nulos no fator de Cholesky nas posições em que na matriz original os elementos são nulos. A observação \ref{obs1} mostra como o preenchimento ocorre.

\begin{remTEMAp}Considere a matriz simétrica definida positiva de ordem $5\times 5$, com o padrão de esparsidade dado à esquerda, em que $\times$ são elementos não nulos. Apenas a parte triangular inferior da matriz está representada. A parte superior pode ser obtida pela transposição da inferior. 
\begin{center}
$\begin{bmatrix}
\times  &        &       &       &     \\
\times  & \times &       &       &      \\
0       & \times & \times &      &      \\
\times  &  0     &    0   & \times &      \\   
0       &    0   &    0   &\times  &\times %   
\end{bmatrix}
\longrightarrow
\begin{bmatrix}
\times &         &       &       &     \\ 
\times  & \times &       &       &      \\
0       & \times & \times &      &      \\
\times  & +      &    +   & \times &      \\   
0       &  0     &    0   &\times  &\times %   
\end{bmatrix}$
\end{center}\label{obs1} 

O cálculo da segunda coluna de $L$ gera um preenchimento na posição $(2,4)$, como ilustrado na matriz à direta, indicado com o sinal $+$ nesta posição. Uma vez calculada a segunda coluna, o cálculo da terceira coluna produzirá um elemento não-nulo na posição $(3,4)$, também indicado pelo símbolo $+$.
\end{remTEMAp}

O preenchimento gerado na fatoração de Cholesky leva à necessidade de mais espaço de armazenamento e maior tempo computacional para a resolução dos sistemas lineares. Logo, seu uso pode se tornar inviável por limitação de tempo e memória, fazendo com que a utilização de métodos iterativos  se torne mais apropriada. Neste contexto, a escolha mais natural de um método para resolver o sistema de equações normais é o método dos gradientes conjugados precondicionado (GCP).

Em geral, a matriz do sistema de equações normais sofre mudanças significativas entre as iterações dos MPI primal-dual e se torna extremamente mal condicionada nas iterações finais, por isso, junto ao método GCP nos MPI é mais adequado a utilização de precondicionadores híbridos compostos por mais de um tipo  de precondicionador, em que cada um é mais apropriado para uma determinada etapa da resolução do problema \cite{GOS13,VOC11}. 

Um precondicionador que se mostra eficiente para as iterações iniciais dos MPI é construído pela fatoração controlada de Cholesky(FCC) \cite{Ca95}, que é um tipo de fatoração incompleta de Cholesky.

Neste tipo de fatoração, ao fatorar uma matriz simétrica definida positiva, é obtida uma matriz triangular inferior $\tilde{L}$, denominada fator incompleto de Cholesky, cujos elementos são calculados como na fatoração de Cholesky, porém  devido a algum tipo de controle de preenchimento, $\tilde{L}$ é mais esparsa que o fator de Cholesky. Na seção \ref{FCC}, a fatoração controlada de Cholesky será descrita mais detatalhadamente.

A consequência de resolver os sistemas lineares que surgem  dos métodos de pontos interiores primal-dual de forma iterativa,  é que a direção de busca em cada iteração é calculada aproximadamente. A direção calculada satisfará um sistema da forma

\begin{equation}\label{sistema} 
\begin{bmatrix}
A & 0 & 0 \\ 
0 & A^{T} & I   \\ 
Z & 0 & X%
\end{bmatrix}%
\begin{bmatrix}
dx \\ 
dy \\ 
dz%
\end{bmatrix}%
=%
\begin{bmatrix}
r_{p} \\ 
r_{d} \\ 
r_{c}%
\end{bmatrix}%
+
\begin{bmatrix}
\tilde{r_1}\\
\tilde{r_2}\\
\tilde{r_3}%
\end{bmatrix}% 
\end{equation}
em que $\tilde{r_1}$, $\tilde{r_2}$ e $\tilde{r_3}$ são componentes de um erro residual.

Se o sistema de equações normais  é resolvido por um método iterativo então a direção $dy$ calculada satisfaz
\begin{center}
$ADA^T dy=  h + \tilde{r}$
\end{center}
em que $h$ é como em (\ref{h}) e  $\tilde{r}$ é um erro residual. Desta forma, o erro residual de (\ref{sistema}) será dado por:
\begin{align}
\tilde{r_1}=\tilde{r}=ADA^Tdy-h\label{r1}\\
\tilde{r_2}=0\label{r2}\\
\tilde{r_3}=0.\label{r3}
\end{align}

Os métodos de pontos interiores que resolvem o sistema de Newton de forma aproximada são conhecidos como métodos de pontos interiores inexatos. Em geral, nestes métodos, uma direção de busca calculada é aceitável se o erro residual satisfaz a condição:
\begin{center}
$\begin{bmatrix}
\tilde{r_{_1}}^k\\
\tilde{r_{_2}}^k\\
\tilde{r_{_3}}^k%
\end{bmatrix}%
\leq\gamma^k\begin{bmatrix}
r^k_{_p} \\ 
r^k_{_d} \\ 
r^k_{_c}%
\end{bmatrix}%
$
\end{center}
em que  $\gamma^k$ é uma sequência uniformemente menor que um.

As análises dos métodos de pontos interiores inexatos mostram que se o erro residual é controlado, ainda é possível calcular uma boa direção de busca em cada iteração \cite{AG09,BS97,Bel98, MON03}.

\newsec{Fatoração Controlada de Cholesky} \label{FCC}

Em uma fatoração incompleta de Cholesky de uma matriz simétrica definida postiva $\mathcal{A}$ de ordem $m$, uma vez calculada uma coluna de $\tilde{L}$ como na fatoração de Cholesky, o controle sobre o preenchimento é feito descartando por algum critério um certo número de elementos não-nulos da coluna. Após este descarte, passa-se para o cálculo da próxima coluna, realizando o mesmo procedimento. Desta forma, é obtida uma matriz triangular inferior $\tilde{L}$, tal que $\mathcal{A}=\tilde{L}\tilde{L}^{T}-R$, em que $R$ é uma matriz denominada matriz resto.


Em geral, as fatorações incompletas de Cholesky são utilizadas como precondicionadores para os métodos de gradientes conjugados precondicionados. Em \cite{DM89} é mostrado que o número de iterações do GCP está quase que diretamente relacionado com a norma da matriz $R$. Isto motivou Campos  \cite{Ca95} a construir a FCC baseada na minimização da norma de Frobenius da matriz $E=\tilde{L}-L$, em que $L$ é o fator de Cholesky da matriz $\mathcal{A}$. De fato, reescrevendo $R$ da forma


\begin{center}
$\begin{array}{rl}
R= & LL^T-\tilde{L}\tilde{L}^T\\
= & LL^T-L\tilde{L}^T+L\tilde{L}^T-\tilde{L}\tilde{L}^T\\
= & L(L-\tilde{L})^T+(L-\tilde{L})\tilde{L}^T,%

\end{array}$
\end{center}
pode-se observar que, $\Vert R \Vert\rightarrow 0$ quando $\Vert E \Vert\rightarrow 0$.


O problema de minimização da norma de Frobenius da matriz $E$ pode ser escrito como
\begin{equation}
min \Vert E \Vert ^2_F= min\left\lbrace\sum_{j=1}^{m} c_j\right\rbrace,
\end{equation}
com $c_j= \sum_{i=1}^{m} \vert l_{ij} - \tilde{l_{ij}} \vert ^2,$ em que $l_{ij}$ e $\tilde{l}_{ij}$, $i,j=1,...,m$ são os elementos de $L$ e $\tilde{L}$, respectivamente. Os elementos $c_j$, por sua vez,  podem ser reescritos como soma de dois somatórios
\begin{equation}
c_j= \sum_{k=1}^{n_j + \eta}\vert l_{i_kj} - \tilde{l}_{i_kj} \vert ^2 + \sum_{k=n_j + \eta +1}^{m}\vert l_{i_kj} \vert ^2,\label{cj}
\end{equation}
em que $n_j$ é o número de elementos não nulos da coluna $j$ de $\mathcal{A}$ e $\eta$ é um parâmetro que determina a quantidade de preenchimento adicional permitida em relação a $n_j$, na coluna $j$ de $\tilde{L}$

A primeira parcela de (\ref{cj}) contém todos os elementos não nulos da coluna $j$ de $\tilde{L}$ e a segunda contém apenas os elementos restantes do fator completo $L$, cujo correspondente em $\tilde{L}$ é nulo. Se consideramos que $\tilde{l}_{ij}\approx l_{ij}$, uma redução da norma de Frobenius de $E$ pode ser obtida por meio da seguinte heurística:

\begin{itemize}
\item Aumentar o valor de $\eta$ permitindo mais preenchimento. Isto faz com que $c_j$, $j=1,..,m$ e consequentemente $\Vert E \Vert _F$ decresça, pois o primeiro somatório em (\ref{cj}) conterá mais elementos e o segundo menos elementos.
\item Para um $\eta $ fixo, escolher os $n_j+\eta$ maiores elementos da coluna $j$ de $L$, para $j=1,...m$, em valor absoluto. Assim, para uma determinada quantidade de armazenamento, um fator incompleto $\tilde{L}$ ótimo é obtido, já que os maiores elementos estarão no primeiro somatório e os menores no segundo.
\end{itemize}

A FCC combina os dois itens acima. Uma vez calculada a coluna $j$ de $\tilde{L}$, ela propõe utilizar um parâmetro $\eta$ tomando o número máximo de elementos não-nulos desta coluna como sendo $k_j=\eta+n_j$. Os $k_j$ elementos maiores em valor absoluto calculados, serão mantidos na coluna $j$ de $\tilde{L}$ e os demais descartados, isto é, serão nulos. %Após este descarte, inicia-se o cálculo da próxima coluna e o procedimento é repetido.

Esta heurística permite prever o armazenamento necessário na fatoração e exercer um controle de preenchimento de acordo com a memória disponível.

Como para um dado problema, a ordem $m$ e o número de elementos não nulos de $\mathcal{A}$ são fixos, o número de elementos não nulos de $\tilde{L}$ depende apenas do parâmetro $\eta$, que por sua vez, pode assumir qualquer valor de $-m$ a $m$. Se $\eta=-m$, a matriz $\tilde{L}$ será diagonal, se $\eta=0$, ela requer o mesmo espaço de armazenamento da matriz $\mathcal{A}$ e se $\eta=m$, a matriz $\tilde{L}$ é o fator de Cholesky completo.

\subsection{Tratamento de pivôs inadequados na fatoração incompleta de Cholesky}\label{TPI}

A existência de uma fatoração incompleta foi estabelecida apenas para certas classes de matrizes, como a classe das $M$-matrizes e a classe das $H$-matrizes, na qual as matrizes diagonalmente dominantes 
estão incluídas \cite{BENZI2002}. Para matrizes simétricas definidas mais gerais, fatorações incompletas de Cholesky podem falhar devido a ocorrência de pivôs nulos ou negativos calculados durante a fatoração. Infelizmente, muitas aplicações importantes lidam com matrizes que não são $M$-matrizes ou $H$-matrizes. 

Pivôs inadequados, são reconhecidos como um problema desde o início do precondicionamento por fatorações incompletas, e inúmeras estratégias foram propostas para remediá-lo.

Uma boa estratégia,  foi sugerida por Manteuffel \cite{Ma80}. Nela, se uma fatoração incompleta de Cholesky de $\mathcal{A}$ falhar devido a pivôs  inadequados, é realizado um incremento global na diagonal de $\mathcal{A}$, antes de se tentar uma nova fatoração. Assim, a fatoração incompleta é aplicada na matriz $\hat{\mathcal{A}}=\mathcal{A}+\alpha\cdot$ diag$(\mathcal{A})$, em que $\alpha >0$ e diag$(\mathcal{A})$ denota a parte diagonal de $\mathcal{A}$. Se a fatoração incompleta falhar novamente, é realizado um incremento em $\alpha$ e a fatoração é reiniciada. O processo é repetido até que uma fatoração incompleta seja calculada com sucesso. 

Esta estratégia é baseada no fato de que existe um valor $\alpha^*$ para o qual a fatoração incompleta de $\hat{\mathcal{A}}$  existe. Pode-se, por exemplo, tomar $\alpha^*$ como sendo o menor $\alpha$ para o qual $\hat{\mathcal{A}}=\mathcal{A}+\alpha\cdot$ diag$(\mathcal{A})$ é diagonalmente dominante. Contudo, bons resultados são geralmente obtidos para valores de $\alpha$ bem menores que $\alpha^*$ e maiores que o menor $\alpha$ para o qual $\hat{\mathcal{A}}$  admite uma fatoração incompleta. A desvantagem desta abordagem é a dificuldade em se determinar um valor de $\alpha$ que não seja muito maior que o necessário para evitar novos reinícios. 

O código da fatoração controlada de Cholesky de Campos utiliza a estratégia de Manteuffel. Bocanegra, em seu precondicionador híbrido \cite{BCO06}, utilizou este código, com $\alpha$ determinado pela fórmula 
\begin{equation}\label{alpha}
\alpha_i =\lambda 2^i,
\end{equation} tomando $\lambda= 5\cdot 10^{-4}$ e $i=1, 2,..., 15$, em que $i$ é o número de reinícios da fatoração.  


\newsec{Nova abordagem na solução direta do sistema de equações normais}\label{NA}

Se por um lado a fatoração de Cholesky do sistema de equações normais oriundos dos MPI primal-dual  gera muitos preenchimentos, por outro, a resolução deste sistema pelo  método dos gradientes conjugados, exige a utilização de precondicionadores sofisticados. Neste trabalho, com o intuito de contornar o problema de preenchimento gerado na fatoração de Cholesky e com a motivação de que uma direção aproximada com boas propriedades pode ser calculada no método de pontos interiores, é proposta uma abordagem que de certa forma combina características dos métodos acima. A abordagem resolve de forma direta sistemas lineares aproximados do sistema de equações normais, mas exercendo um certo controle sobre o preenchimento. Isto é feito substituindo a fatoração de Cholesky da matriz do sistema de equações normais, pela fatoração controlada de Cholesky.

Desta forma, ao invés de resolver em cada iteração o sistema linear $LL^Tdy=h$, onde $L$ é o fator de Cholesky da matriz $ADA^{T}$, é resolvido ainda diretamente, um sistema linear da forma
\begin{equation}
\tilde{L}\tilde{L}^T\tilde{dy}=h,\label{sol_aprox}
\end{equation}  
em que $\tilde{L}$ é o fator incompleto da matriz $ADA^{T}$, calculado por meio da fatoração controlada de Cholesky. Assim, a cada iteração do método é encontrada uma solução $\tilde{dy}$  que é uma aproximação da solução original, isto é, $\tilde{dy}$ satisfaz o sistema

\begin{equation}
ADA^T\tilde{dy}=h+\tilde{r},
\end{equation}
em que $\tilde{r}=ADA^T\tilde{dy}-\tilde{L}\tilde{L}^T\tilde{dy}$.

Se $\tilde{r}$ é controlado, o método proposto é de certa forma um método de pontos interiores inexato, no sentido de que uma direção aproximada é calculada. Entretanto, difere dos métodos inexatos usuais que resolvem o sistema de equações normais. Os métodos inexatos usuais utilizam um método iterativo para resolver o sistema de equações normais ``exato'', já a abordagem proposta resolve de forma direta um sistema aproximado do sistema de equações normais. Assim, enquanto a cada iteração dos métodos inexatos usuais, muitas iterações do método iterativo são feitas, até que o resíduo (\ref{r1}) satisfaça a condição de convergência, o método proposto realiza uma única vez substituição regressiva e progressiva para resolver os sistemas lineares com as matrizes triangulares $\tilde{L}$ e $\tilde{L}^T$.

O sucesso de um MPI inexato depende de tornar o resíduo que surge do lado direito da equação de Newton menor do que uma tolerância aceitável para atingir a convergência. 

Observamos que, tomando a norma dois de $\tilde{r}$, obtém-se
\begin{center}
$\begin{array}{rl}
\Vert \tilde{r}\Vert _2 & =\Vert ADA^T\tilde{dy}-h\Vert _2\\
 &= \Vert LL^T\tilde{dy}-\tilde{L}\tilde{L}^T\tilde{dy}\Vert _2 \\
 &\leq\Vert LL^T-\tilde{L}\tilde{L}^T\Vert _F \Vert\tilde{dy}\Vert _2\\
 &=\Vert L(L-\tilde{L})^T+(L-\tilde{L})\tilde{L}^T\Vert _F \Vert\tilde{dy}\Vert _2\\
 &\leq(\Vert L\Vert _F \Vert(L-\tilde{L})^T\Vert _F +\Vert L-\tilde{L}\Vert _F \Vert\tilde{L}^T\Vert _F )\Vert\tilde{dy}\Vert _2\\
 & =\Vert L-\tilde{L}\Vert _F(\Vert L\Vert _F + \Vert\tilde{L}^T\Vert _F) \Vert\tilde{dy}\Vert _2.
 \end{array}$
\end{center}
Então, se $\Vert L\Vert _F$, $ \Vert\tilde{L}^T\Vert _F$ e $\Vert\tilde{dy}\Vert _2$ forem limitadas, temos que $\Vert \tilde{r} \Vert_2 \rightarrow 0$, quando $\Vert L-\tilde{L} \Vert _F \rightarrow 0$. Logo, na abordagem proposta é desejável minimizar a norma de Frobenius de $E=L-\tilde{L}$, o que justifica a escolha da FCC para a construção do fator incompleto.

Além disso, o parâmetro de preenchimento $\eta$ da FCC pode ser controlado convenientemente, de forma a iniciar com um valor que permita na primeira iteração, uma matriz $\tilde{L}$ o mais esparsa possível. À medida que o MPI avança nas iterações, o valor de $\eta$ aumenta progressivamente, permitindo um aumento gradual da densidade de $\tilde{L}$ ao longo das iterações, até atingir a fatoração de Cholesky completa. Ao mesmo tempo espera-se que a convergência do método seja alcançada. Desta forma, a resolução do sistema se torna mais rápida nas iterações iniciais, já que haverá uma economia de tempo tanto no cálculo dos elementos do fator incompleto, quanto na resolução dos dois sistemas triangulares.

\newsec{Experimentos Computacionais}

A abordagem proposta foi incorporada ao código PCx \cite{JSMS99} que é uma implementação do método primal-dual preditor-corretor. Todas as rotinas deste código são implementadas em Linguagem C, exceto a parte responsável pela solução dos sistemas lineares que é feita aplicando a fatoração de Cholesky no sistema de equações normais, utilizando o código da fatoração de Cholesky esparsa projetada por Ng e Peyton \cite{NP93}, em FORTRAN77. Na implementação do código com a nova abordagem, que será denominada código PCx modificado (PCx$_{mod}$), foi admitida a existência de duas fases. Na primeira fase são calculados fatores de Cholesky incompletos, utilizando um código da FCC implementado em FORTRAN77 que foi fornecido pelo seu autor \cite{Ca95}, mas com algumas adaptações. Na segunda fase, são calculadas fatorações de Cholesky completas utilizando o código de Ng e Pyton. Isto porque, o código da FCC não é eficiente para o cálculo do fator completo, pois necessita de operações adicionais e não explora o fato de que a estrutura esparsada da fatoração final é conhecida.

A mudança para a segunda fase é feita quando o número de elementos não nulos do fator incompleto calculado é $95\%$ do número de elementos não nulos do fator de Cholesky completo\footnote{O número de elementos não nulos do fator de Cholesky completo é calculado pelo código de Ng e Peyton na fase incial, antes de iniciar as iterações.} ou quando não está sendo obtido  progresso com a fatoração controlada de Cholesky. O progresso de uma iteração para outra é medido através do quociente $\rho=\dfrac{\mu^{k-1}}{\mu^{k}}$, em que $\mu^{k}=(x^k)^tz^k/n $ e $\mu^{k-1}$ são os \textit{gap} de dualidade da iteração atual e anterior, respectivamente. Na implementação foi considerado que se $\rho\geq 0,99$, nenhum progresso está sendo obtido, assim a mudança para a segunda fase é feita na iteração imediatamente seguinte.

Ainda na primeira fase, uma vez determinado um valor de $\eta$ inicial, a partir da segunda iteração, seu valor é atualizado ou não, em função de $\rho$, da seguinte forma: 
\begin{center}
$\eta=\left\lbrace
\begin{array}{ll}
\eta, & \text{se } \rho < 0,3\\
\eta+a\rho, & \text{se } 0,3 \leq \rho \leq 0,7\\
\eta+b\rho, & \text{se } \rho \geq 0,7%
\end{array}\right. $
\end{center}\label{atualizacao_eta}
em $a$ e $b$ são constantes, com $a<b$.

Desta forma, se um bom progresso é obtido de uma iteração para outra, isto é, se há uma redução significativa do \textit{gap} de dualidade, o valor do parâmetro $\eta$ não muda de uma iteração para outra. E se pouco progresso é obtido na redução do \textit{gap} de dualidade, é realizado um incremento maior em $\eta$.

Na seção \ref{TPI}, foi mencionado o problema do surgimento de pivôs inadequados nas fatorações incompletas de Cholesky. No método proposto, para contornar este problema, foi utilizada a mesma técnica que o código PCx utiliza quando pivôs inadequados surgem na fatoração de Cholesky. Um pivô $\tilde{L_{ii}}$ menor que $10^{-8}$ é substituído por $10^{128}$ e a fatoração simplesmente continua a partir daí, evitando reiniciar a fatoração. Isso faz com que os elementos da coluna $i$ de $\tilde{L}$ sejam essencialmente nulos, replicando técnicas utilizadas para fatoração de matrizes semidefinidas positivas \cite{Hig96}.

Os testes apresentados nesta seção foram realizados com os problemas de fluxo multiproduto do tipo PDS da biblioteca Kennington \cite{KENN90}. A tabela \ref{tabprob1} exibe em suas colunas, nesta ordem, o nome dos problemas, o número de linhas, o número de colunas, a quantidade de elementos não nulos da matriz $AA^T$ ($|AA^T|$) e  a densidade do fator de Cholesky desta última matriz.
\vspace{-0.5cm}

\begin{table}[H]
\centering
\caption{Problemas testes}
\vspace{0.2cm}
\begin{tabular}{lrrrc}
\hline
\hline
Problema & Linhas & Colunas & $|AA^T|$ & Densidade Chol
\\
\hline 
pds-10 & 16558 & 49932 & 77866 & 0,014 
\\ 
pds-20 & 33874 & 108175 & 170973 & 0,130 
\\
pds-30 & 49944  & 158489  &251580 & 0,013 
\\
pds-40 &	 64265 &	214385 & 342483   &0,013 
\\
pds-50 & 83060 &  275814  & 432371& 0,013 
\\
pds-60 &	 96503 & 332862 & 524563 & 0,012 
\\
pds-70 &  114944 & 390005 & 607375& 0,013 
\\
pds-80 &	 126120 & 430800 & 677160 &  0,012 
 \\
pds-90 &  142823 & 475448 & 741153& 0,011 
\\ 
pds-100 & 156243 & 514577 & 788421 & 0,010 
\\
\hline
\hline 
\end{tabular} \label{tabprob1}
\end{table}

Nos testes foram utilizados os valores $0$, para o valor inicial de $\eta$, $a=10$ e $b=25$, sendo estes os melhores valores obtidos para os problemas PDS, dentre alguns valores testados.

A tabela \ref{resultados}, exibe uma comparação, na segunda e na terceira coluna, entre o número de iterações dos códigos PCx e PCx$_{mod}$, respectivamente. Na quarta coluna temos o número de iterações da primeira fase do código, ou seja, o número de iterações em que se utilizou um fator incompleto para calcular uma direção de busca. Nas duas últimas colunas da tabela, pode-se comparar o tempo de processamento, em segundos, para a resolução dos problemas pelos dois códigos. A última linha da tabela exibe a soma total de cada coluna.

\begin{table}[H] 
\scalefont{0.91}
\centering
\caption{Comparação do número de iterações e do tempo de processamento}
\vspace{0.2cm}
\begin{tabular}{lrrrrrrcrrrr}
 
\hline \hline & \multicolumn{3}{c}{Iteração} & &\multicolumn {2}{c}{ Tempo} \\ 
  
\cline{2-4} \cline{6-7} \raisebox{0.5ex}{ Problema }& PCx & PCx$_{mod}$ & FCC && PCx & PCx$_{mod}$ \\ 
  
\hline 

pds-10 & \textbf{40} &  55   & 20  && 30,61 & \textbf{17,79}\\
pds-20 & \textbf{52} &  66   & 29  && 177,01 & \textbf{140,09}\\
pds-30 & \textbf{62} &  66   & 26  && 563,14 & \textbf{338,26}\\%ok
pds-40 & 69 & \textbf{67}& 25 &&1826,69&\textbf{978,52} \\ 
pds-50 & \textbf{68} & 72    & 30 && 2760,82  & \textbf{1700,20}\\
pds-60 & \textbf{69} & 77 & 34 && 4495,51 & \textbf{2883,73} \\ %ok
pds-70 & \textbf{72} & 82 & 39  &&7425,77  & \textbf{4562,90}\\
pds-80 & \textbf{74} & 75 & 33 && 10631,58& \textbf{6083,00} \\ %ok
pds-90 & \textbf{75} &  83   & 45  &&  13348,96 & \textbf{6845,25} \\
pds-100 & \textbf{80} &  84   &  11 && 15724,24 & \textbf{14222,93}\\
\hline 
Total &  661  & 727  &  292  & & 56984,33   & 37772,67
 \\
\hline \hline \end{tabular} \label{resultados} 
\end{table}

Com exceção do problema \textit{pds-40}, houve um aumento no número de iterações dos demais problemas, somando ao todo uma aumento de 66 iterações, como pode ser visto na tabela \ref{resultados}. Entretanto, em contrapartida, houve uma redução significativa no tempo de processamento de todos os problemas, totalizando uma redução de 19211,66 segundos.

A FCC foi utilizada em $40\%$ do total de iterações e a redução total do tempo de processamento foi de $33,7\%$. A grande redução obtida no tempo se deveu à resolução dos sistemas lineares mais esparsos nas iterações iniciais. 

A tabela \ref{resultados1} permite analisar, para cada problema, o quanto o fator incompleto de Cholesky calculado na primeira iteração, é menos denso que o fator de Cholesky. Isto pode ser observado, comparando o número de elementos não nulos do fator de Cholesky $L$ e do fator incompleto da primeira iteração, denotado por  $\tilde{L_0}$, exibidos na segunda e terceira coluna da tabela, respectivamente. Note que neste caso como $\eta$ inicial é igual a $0$, o número de elementos não nulos da matriz  $\tilde{L_0}$ é igual ao número de elementos não nulos da matriz $AA^T$. Os valores do quociente 
\begin{center}
$\kappa=\dfrac{|\tilde{L_0}|}{|L|}$,
\end{center}
exibidos na quarta coluna, nos fornecem uma melhor visão da esparsidade de $\tilde{L_0}$ em relação a $L$, pois, se $\kappa$ é próximo de zero, então $\tilde{L_0}$  é muito mais esparso que $L$ e quando $\kappa$ é próximo de um, $\tilde{L_0}$ é quase o fator de Cholesky. 
\begin{table}[H]
\scalefont{0.95}
\centering
\caption{Comparação entre Cholesky Completo e o Incompleto}
\vspace{0.2cm}

\begin{tabular}{lrrc}
\hline
\hline
Problema & $|L|$ & $|\tilde{L_0}|$ &$\kappa$ 
\\
\hline


pds-10 & 1687660 & 79866  & 0,047 
\\
pds-20 & 7089645 & 170973  & 0,024 
\\
pds-30 & 14643812 & 251580 & 0,017 
\\
pds-40 &	 28195225 &	342483 & 0,012 
\\
pds-50 & 41767191 & 432371  & 0,010 
\\
pds-60 &	 58118583 & 5224563& 0,009 
\\
pds-70 & 78724946 & 607375 & 0,007 
\\
pds-80 &	 94270275 & 677160  & 0,008 
 \\ 
pds-90 & 110971961 & 741153 & 0,007 
\\
pds-100 &120240013 & 788421 & 0,006 
\\
\hline
\hline 
\end{tabular} \label{resultados1} 

\end{table}

Pode-se observar na tabela \ref{resultados1} que o primeiro fator incompleto das matrizes dos sistemas de equações normais dos problemas do tipo PDS são extremamente esparsas em comparação ao fator de Cholesky completo, o que justifica o redução no tempo de processamento total apesar do aumento do número de iterações.

\subsection{Erro residual}

Na prática, os métodos de pontos interiores inexatos usuais, que resolvem o sistema de equações normais, utilizam um método iterativo para resolver este sistema. O critério de parada do método iterativo é baseado no controle  do erro residual, assim, as iterações são realizadas, até que o erro residual atinja uma condição pré-estabelecida para que a convergência ocorra.

Na abordagem proposta, em uma iteração de métodos de pontos interiores, não é possível a priori, controlar o erro residual, podendo ocorrer de ele ser maior do que a condição pré-estabelecida para a convergência. Entretanto, dado que erro residual $\tilde{r}$ no lado direito do sistema de equações normais satisfaz a condição

\begin{center}
$\Vert \tilde{r}\Vert _2 \leq \Vert L-\tilde{L}\Vert _F(\Vert L\Vert _F + \Vert\tilde{L}^T\Vert _F) \Vert\tilde{dy}\Vert _2$,
\end{center}
caso seja necessário, pode-se aumentar a densidade da matriz $\tilde{L}$ de forma que ela se aproxime mais da matriz $L$, até que $\Vert \tilde{r}\Vert _2$ satisfaça a condição de convergência. Na pior das hipóteses, pode ser tomado $\tilde{L}=L$.

Assim, calculada uma direção $\tilde{dy}$ aproximada em uma iteração de método de pontos interiores, é verificado se ela satisfaz a condição de convergência. Se isto não ocorre, descarta-se a direção calculada, calcula-se uma nova matriz $\tilde{L}$ permitindo mais preenchimentos e uma nova direção $\tilde{dy}$ é calculada utilizando esta matriz.

Em \cite{Ko00}, Korzak propõe um método de pontos interiores primal-dual infactível e inexato baseado no método proposto por Kojima, Mizuno e Yoshise em \cite{KMY89}. No método de Korzak a convergência é garantida se as componentes do erro residual no sistema (\ref{sistema})  satisfazem as condições:

\begin{align}
\Vert \tilde{r_{_1}}^k\Vert \leq (1-\tau _1)\Vert r^k_{_p}\Vert;\label{cond1}\\
\Vert \tilde{r_{_2}}^k\Vert \leq (1-\tau _2)\Vert r^k_{_d}\Vert;\label{cond2}\\
\Vert \tilde{r_{_3}}^k\Vert \leq (1-\tau _3)\mu ^k,\label{cond3} 
\end{align}

com $\tau_1\in (0,1]$,  $\tau_2\in (0,1]$ e  $\tau_3\in (0,1]$. 
Na prática, é utilizado $\tau_1=0,95$,  $\tau_2=0,95$ e  $\tau_3=0,049$.

Baseado no trabalho de Korzak, foram realizados alguns testes computacionais para verificar o comportamento do erro residual na abordagem proposta.  Neste caso, as condições (\ref{cond2}) e (\ref{cond3})  são atendidas, já que $\tilde{r_2}=0$ e  $\tilde{r_3}=0$. Assim, como $\tilde{r}=\tilde{r_1}$, é verificado apenas o valor da norma  $\left\Vert\frac{\tilde{r_{_1}}^k}{r^k_{_p}}\right\Vert$ no processamento de alguns problemas, nas iterações em que a fatoração incompleta de Cholesky foi realizada. A tabela \ref{norma} exibe estes valores para as quinze primeiras iterações de cinco dos problemas testes.

\begin{table}\centering
\caption{ Valor da norma $\left\Vert\frac{\tilde{r_{_1}}^k}{r^k_{_p}}\right\Vert$ nas quinze primeiras iterações}
\begin{tabular}{c|l|l|l|l|l} 
\hline\hline

Iteration & pds-10 & pds-30 & pds-50 & pds-70 & pds-90  \\
\hline
1 &$4\cdot 10^{-2}$ &$4\cdot 10^{-2}$ &$3\cdot 10^{-2}$ & $3\cdot 10^{-2}$&$2\cdot 10^{-2}$\\
2 &$1\cdot 10^{-4}$ &$6\cdot 10^{-3}$ &$1\cdot 10^{-3}$ & $4\cdot 10^{-3}$&$5\cdot 10^{-3}$ \\
3 &$3\cdot 10^{-4}$ &$3\cdot 10^{-4}$ &$5\cdot 10^{-4}$ & $1\cdot 10^{-4}$ &$2\cdot 10^{-4}$ \\
4 & $2\cdot 10^{-5}$&$6\cdot 10^{-5}$ &$2\cdot 10^{-5}$ & $3\cdot 10^{-5}$&$2\cdot 10^{-5}$\\
5 &$3\cdot 10^{-4}$ &$1\cdot 10^{-4}$ &$1\cdot 10^{-4}$ & $1\cdot 10^{-5}$& $2\cdot 10^{-5}$ \\
6 & $3\cdot 10^{-4}$& $2\cdot 10^{-4}$&$1\cdot 10^{-5}$ & $1\cdot 10^{-5}$&$5\cdot 10^{-5}$ \\
7 &$9\cdot 10^{-4}$ &$3\cdot 10^{-4}$ & $1\cdot 10^{-4}$& $2\cdot 10^{-5}$&$1\cdot 10^{-4}$ \\
8 & $2\cdot 10^{-4}$& $4\cdot 10^{-4}$& $2\cdot 10^{-4}$& $1\cdot 10^{-4}$&$1\cdot 10^{-4}$ \\
9 &$1\cdot 10^{-4}$ & $4\cdot 10^{-4}$& $5\cdot 10^{-4}$& $1\cdot 10^{-3}$&$5\cdot 10^{-4}$ \\
10 &$1\cdot 10^{-4}$ &$4\cdot 10^{-4}$ &$1\cdot 10^{-3}$& $1\cdot 10^{-4}$&$6\cdot 10^{-4}$ \\
11 &$1\cdot 10^{-4}$ & $1\cdot 10^{-4}$& $9\cdot 10^{-4}$& $7\cdot 10^{-4}$&$7\cdot 10^{-4}$ \\
12 & $1\cdot 10^{-3}$&$1\cdot 10^{-4}$ & $1\cdot 10^{-3}$& $5\cdot 10^{-4}$&$6\cdot 10^{-4}$ \\
13 &-                 &$1\cdot 10^{-4}$ & $1\cdot 10^{-3}$& $5\cdot 10^{-4}$&$4\cdot 10^{-3}$ \\
14 & -                & $1\cdot 10^{-4}$&$1\cdot 10^{-4}$ &$9\cdot 10^{-4}$&$1\cdot 10^{-3}$ \\
15 & -                & $1\cdot 10^{-3}$& $8\cdot 10^{-4}$&$1\cdot 10^{-3}$&$1\cdot 10^{-3}$\\ 

\hline\hline 
\end{tabular}  \label{norma}
\end{table}
Na tabela pode ser verificado que  o valor de  $\left\Vert\frac{\tilde{r_{_1}}^k}{r^k_{_p}}\right\Vert$ é menor do que $0,05$ em todas as iterações. Assim a condição de convergência(\ref{cond1}) é atendida nestas iterações.

A desvantagem da abordagem proposta é que uma direção aproximada calculada pode não satisfazer a condição de convergência, entretanto nos testes computacionais verificamos que nas quinze primeiras iterações dos problemas testados, nenhuma direção inadequada foi calculada. Mas caso isso ocorra, a direção inadequada é descartada e aumentando a densidade da matriz calculada pela fatoração incompleta, é possível calcular uma nova direção mais adequada. 

\newsec{Conclusões}


Neste trabalho foi apresentado uma nova abordagem para a solução do sistema de equações normais oriundos dos problemas de fluxo multiproduto em métodos de pontos interiores.

O método proposto é de certa forma um método de pontos interiores inexato, no sentido de que uma direção aproximada é calculada em cada iteração. Mas diferente dos métodos inexatos usuais que resolvem o sistema de equações normais ``exato'' por meio de um método iterativo, a abordagem proposta resolve de forma direta um sistema aproximado do sistema de equações normais exercendo um controle sobre o preenchimento.

Os testes computacionais, mostraram que a abordagem proposta é bastante eficiente para a resolução dos problemas de fluxo multiproduto do tipo PDS, sendo competitivo com o código PCx. 

A desvantagem do método é que não é possível prever a priori, se uma direção calculada satisfaz as condições de convergência, entretanto esta direção pode ser descartada e uma uma nova direção, mais adequada, pode ser eventualmente calculada.

Como trabalho futuro iremos extender este método para outros tipos de problemas de programação linear.

%
%\newsec{Agradecimentos}
%
%Os autores agradecem à FAPESP e ao CNPq pelo suporte financeiro deste trabalho.
%
%

\begin{abstract}
{\bf Abstract}. One of the approaches used to solve the linear system that arises at each iteration in the primal-dual interior points methods is  reduce it to an equivalent symmetrical positive defined linear system, known as normal equation system  and apply the  factorization in the system matrix.
The major disavantage of this approach is the filling generated during factorization, which may turn their use impractical for limited time and memory. In order, to deal with the fill-in problem in the Cholesky factorization,  in this paper, we proposed an approach that solve directly approximate  normal equations linear systems  derived  from the multiproduct commodity problems and apply some control over the filling.
\end{abstract}


\begin{thebibliography}{8}

\bibitem{AG09} G. Al-Jeiroudi and J. Gondzio, About the
``Convergence Analysis of the Inexact Infeasible Interior-Point Method for Linear Optimization'',
Journal of Optimization Theory and Applications, {\bf 141} (2009), pp. 231--247.

\bibitem{BS97} V. Baryamureeba and T. Steihaug, ``On the convergence of an inexact
primal-dual interior point method for linear programming'', in Lecture Notes in
Computer Science, Springer, Berlin/Heidelberg, 2006.

\bibitem{Bel98} S. Bellavia, ``An inexact interior point method'',
Journal of Optimization Theory
and Applications, {\bf 96} (1988), pp. 199--121.

\bibitem{BENZI2002} M. Benzi, ``Preconditioning techniques for large linear systems: A survey, Journal
of Computational Physics'', {\bf 182}, No 2 (2002), pp. 418--477.

\bibitem{BCO06} S. Bocanegra, F. F. Campos, and A. R. L. Oliveira, ``Using a hybrid
preconditioner for solving large-scale linear systems arising from interior point
methods'', Computational Optimization and Applications, {\bf 36} (2007), pp. 149--164. Special issue on ``Linear Algebra Issues arising in Interior Point Methods''.

\bibitem{Ca95} F. F. Campos, ``Analysis of Conjugate Gradients - type methods for solving linear
equations'', PhD thesis, Oxford University Computing Laboratory, Oxford,
1995.

\bibitem{JSMS99} J. Czyzyk, S. Mehrotra, S. J. Wright, and M. Wagner, ``PCx: An
interior-point code for linear programming'', Optimization Methods and Software,
{\bf 11}, No 1 (1999), pp. 397--430.

\bibitem{DM89} I. S. Duff and G. A. Meurant, ``The effect of ordering on preconditioned
conjugate gradients'', BIT, 29 (1989), pp. 635--657.

\bibitem{GOS13} C. T. L. S. Ghidini, A. R. L. Oliveira, and D. C. Sorensen, ``Computing a
hybrid preconditioner approach to solve the linear systems arising from interior
point methods for linear programming using the gradient conjugate method'',
Annals of Management Science, (2013).

\bibitem{Hig96} N. Highan, ``Acuraccy and Stability of Numerical Algorithms'', SIAM Publications,
Philadelphia, 1996.

\bibitem{KENN90} J. L. Kennington, S. Niemi, and S. J. Wichmann, ``An empirical evaluation
of the korbx algorithms for military airlift applications'', Operations Research,
{\bf 28}, pp. 240--248.

\bibitem{KMY89} M. Kojima, S. Mizuno, and A. Yoshise, ``A primal-dual interiorpoint
method for linear programming. progress in mathematical programming'',
interior-point and related methods, scSpring-Verlag, New York, (1989), pp. 29--
47.

\bibitem{Ko00} J. Korzak, Convergence analysis of inexact infeasible-interior-point algorithms
for solving linear programming problems'', Siam J. Optim, {\bf 11} (2000),
pp. 133--148.

\bibitem{Ma80} T. A. Manteuffel, ``An incomplete factorization technique for positive definite
linear systems'', Math. Comp., {\bf 34} (1980), pp. 473--497.

\bibitem{MON03} R. D. S. Monteiro and J. W. O'Neal, ``Convergence analysis of long-step
primal-dual infeasible interior point LP algorithm based on iterative linear solvers'',
Georgia Institute of Technology, 2003.

\bibitem{NP93} E. Ng and B. P. Peyton, ``Block sparse Cholesky algorithms on advanced uniprocessors
computers'', SIAM Journal on Scientific Stat. Computing, {\bf 14} (1993),
pp. 1034--1056.

\bibitem{Va96} R. J. Vanderbei, ``Linear Programming - Foundations and Extensions'', Kluwer
Academics Publishers, Boston, USA, 1996.

\bibitem{VOC11} M. I. Velazco, A. R. L. Oliveira, and F. F. Campos, ``Heuristics for
implementation of a hybrid preconditioner for interior-point methods'', Pesquisa
Operacional, {\bf 34} (2011), pp. 2553--2561.

\bibitem{WR96} S. J. Wright, ``Primal-Dual Interior-Point Methods'', SIAM Publications,
SIAM, Philadelphia, PA, USA, 1996.


\end{thebibliography}


\end{document}
\newpage
$ \  \  $  \thispagestyle{myheadings}  \markboth{      }{   }
