Home > source > @penlab > solve_chol.m

# solve_chol

## PURPOSE

(Modified) Newton method via Cholesky factorization (with inertia control)

## SYNOPSIS

function [dir,nFlag]=solve_chol(obj,matrix,rhs)

## DESCRIPTION

     (Modified) Newton method via Cholesky factorization (with inertia control)
     Symmetric approximate minimum degree permutation is possible to employ
     Solves the system "matrix*dir=rhs"
     OUTPUT
       nFlag ... 0..OK, 1..Cholesky/inertia failed
     'obj' is used only for option settings (and to store last perturbation)

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- penlab PenLab (Pennon Laboratory for Matlab, previously PennonM)

## SOURCE CODE

    0001 function [dir,nFlag]=solve_chol(obj,matrix,rhs)
    0002 % (Modified) Newton method via Cholesky factorization (with inertia control)
    0003 % Symmetric approximate minimum degree permutation is possible to employ
    0004 % Solves the system "matrix*dir=rhs"
    0005 % OUTPUT
    0006 %   nFlag ... 0..OK, 1..Cholesky/inertia failed
    0007 % 'obj' is used only for option settings (and to store last perturbation)
    0008
    0009 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0010
    0011
    0012   start_time=cputime;
    0013   [n m] = size(matrix);
    0014   mat_nnz=nnz(matrix);
    0015
    0016   LMUPDATE = obj.allopts.pert_update;
    0017   LMLOW = obj.allopts.pert_min;
    0018   PERT_TRY_MAX = obj.allopts.pert_try_max;
    0019   ORDERING = obj.allopts.chol_ordering;
    0020   %LMFASTER = obj.allopts.pert_faster; % && issparse(matrix); ??
    0021   LMFASTER = 0; % && issparse(matrix); ??
    0022
    0023   if (ORDERING && issparse(matrix))
    0024     if (~isempty(obj.chol_factperm))
    0025       obj.print(4,Inf,'Reusing ordering');
    0026       p=obj.chol_factperm;
    0027     else
    0028       obj.print(4,Inf,'Computing ordering');
    0029       p=symamd(matrix);
    0030       obj.chol_factperm=p;
    0031     end
    0032   else
    0033     p=1:n;
    0034   end
    0035   rev(p)=1:n;  % inverse permutation
    0036
    0037   matrix=matrix(p,p);
    0038   [factor, status] = chol(matrix);
    0039   chol_try=0;
    0040   lambda=max(obj.chol_lmlast/LMUPDATE, LMLOW);  % lambda to try if it fails now
    0041   lambda_last=0.;
    0042   while (status~=0 && chol_try<PERT_TRY_MAX)
    0043     obj.print(4,Inf,'Chol fact failed (%d), new pert: %e',status,lambda);
    0044     if (LMFASTER)
    0045       % find negative curvature vector from the last factorization
    0046       %k = status-1;
    0047       %R11 = factor(1:status-1,1:status-1);
    0048       %R12 = factor(1:status-1,status:n);
    0049       %v_nc = [R11 \ R12(:,1); -1; zeros(n-k-1,1)];
    0050       R12e1 = factor(1:status-1,status);
    0051       v_nc = [factor(1:status-1,1:status-1) \ R12e1; -1; zeros(n-status,1)];
    0052       v_nc = v_nc/norm(v_nc);      % negative curvature vector
    0053       % perturbated matrix must be positive definite, e.i., v_nc'*(matrix+lambda*I)*v_nc>0
    0054       % it gives: lambda > -v_nc'*matrix*v_nc since v_nc is normed
    0055       negcrv = -v_nc'*matrix*v_nc;
    0056       cnt=0;
    0057       while (lambda-lambda_last <= negcrv)
    0058         lambda=lambda*LMUPDATE;
    0059         cnt=cnt+1;
    0060       end
    0061       obj.print(5,Inf,'Chol: pert_faster: skipped %d steps: pert=%e -> %e',cnt,lambda_last,lambda);
    0062     end
    0063     matrix=matrix + (lambda-lambda_last)*speye(n,n);
    0064     [factor, status] = chol(matrix);
    0065     chol_try = chol_try+1;
    0066     lambda_last=lambda;
    0067     lambda=lambda*LMUPDATE;
    0068   end
    0069
    0070   if (chol_try>0)
    0071     obj.chol_lmlast=lambda_last;
    0072   end
    0073
    0074   % update statistics
    0075   time_fact = cputime - start_time;
    0076   obj.stats_time_fact_last=obj.stats_time_fact_last+time_fact;
    0077   obj.initer_last = obj.initer_last+chol_try+1;
    0078
    0079   if (status~=0)
    0080     obj.print(3,Inf,'Chol fact failure (%d), last pert: %e, giving up',status,lambda_last);
    0081     nFlag = 1;
    0082     dir=zeros(n,1);
    0083     return;
    0084   end
    0085
    0086   dir=factor \ (factor' \ rhs(p));
    0087   dir=dir(rev);
    0088   time_total=cputime - start_time;
    0089   nFlag=0;
    0090
    0091   obj.print(4,Inf,'Chol fact OK in %fs, total %fs, no pert=%i, pert=%.4e, nnz=%d (dim %dx%d)',time_fact, time_total, chol_try, lambda_last, nnz(factor),n,n);
    0092   if (lambda_last>0)
    0093     obj.print(3,4,'Chol fact: OK, pert=%.4e',lambda_last);
    0094   end
    0095
    0096