%% 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