Home > source > @penlab > solvekkt_ldl.m

# solvekkt_ldl

## PURPOSE

A direct solver based on LDL() function in Matlab (version 7.6, R2008a++)

## SYNOPSIS

function [x_dir,ueq_dir,nFlag,pert]=solvekkt_ldl(obj,H,A,rhs1,rhs2)

## DESCRIPTION

     A direct solver based on LDL() function in Matlab (version 7.6, R2008a++)
     for indefinite KKT system [H,A; A',0]
     with inertia control and diagonal perturbations if needed
     INPUT
       H, A  ... matrices nxn, nxm respectively forming the system K=[H,A;A',0]
       rhs1,2... vectors nx1, mx1, right hand side
     OUTPUT
       x_dir, ueq_dir ... primal & dual reached solution
       nFlag ... flag (0..OK, 1..inertia control/ldl() failed)
       pert ... optional parameter - used perturbation of (1,1) block to make
          the kkt matrix inertia right

     'obj' is used for option settings

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

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

## SUBFUNCTIONS

- function [e_pos, e_neg, e_zero] = inertia(D)

## SOURCE CODE

    0001 function [x_dir,ueq_dir,nFlag,pert]=solvekkt_ldl(obj,H,A,rhs1,rhs2)
    0002 % A direct solver based on LDL() function in Matlab (version 7.6, R2008a++)
    0003 % for indefinite KKT system [H,A; A',0]
    0004 % with inertia control and diagonal perturbations if needed
    0005 % INPUT
    0006 %   H, A  ... matrices nxn, nxm respectively forming the system K=[H,A;A',0]
    0007 %   rhs1,2... vectors nx1, mx1, right hand side
    0008 % OUTPUT
    0009 %   x_dir, ueq_dir ... primal & dual reached solution
    0010 %   nFlag ... flag (0..OK, 1..inertia control/ldl() failed)
    0011 %   pert ... optional parameter - used perturbation of (1,1) block to make
    0012 %      the kkt matrix inertia right
    0013 %
    0014 % 'obj' is used for option settings
    0015
    0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0017
    0018   start_time = cputime;
    0019   [n m] = size(A);
    0020
    0021   LMUPDATE = obj.allopts.pert_update;
    0022   LMLOW = obj.allopts.pert_min;
    0023   PERT_TRY_MAX = obj.allopts.pert_try_max;
    0024   pivot = obj.allopts.ldl_pivot;
    0025
    0026   % LDL() with inertia control
    0027   block22=sparse(m,m);
    0028
    0029   [L,D,perm] = ldl([H,A;A',block22],pivot,'vector');
    0030   [e_pos, e_neg, e_zero] = inertia(D);
    0031
    0032   pert_try=0;
    0033   lambda=max(obj.chol_lmlast/LMUPDATE, LMLOW);  % lambda to try if it fails now
    0034   lambda_last=0.;
    0035
    0036   while ((e_zero>0 || e_neg>m) && pert_try<PERT_TRY_MAX)
    0037     obj.print(4,Inf,'Wrong inertia: %d,%d,%d, new pert: %e',e_pos,e_neg,e_zero,lambda);
    0038     H = H + (lambda-lambda_last)*speye(n,n);
    0039     [L,D,perm] = ldl([H,A;A',block22],pivot,'vector');
    0040     [e_pos, e_neg, e_zero] = inertia(D);
    0041
    0042     pert_try = pert_try+1;
    0043     lambda_last=lambda;
    0044     lambda=lambda*LMUPDATE;
    0045   end
    0046
    0047   if (pert_try>0)
    0048     obj.chol_lmlast=lambda_last;
    0049   end
    0050
    0051   % update statistics
    0052   time_factor = cputime-start_time;
    0053   obj.stats_time_fact_last=obj.stats_time_fact_last+time_factor;
    0054   obj.initer_last = obj.initer_last+pert_try+1;
    0055   if (nargout>=4)
    0056     pert=lambda_last;
    0057   end
    0058
    0059   if (e_zero~=0 || e_neg~=m)
    0060     obj.print(3,Inf,'Inertia control failure, wrong inertia: %d,%d,%d, last pert: %e',e_pos,e_neg,e_zero,lambda_last);
    0061     nFlag = 1;
    0062     x_dir=zeros(n,1);
    0063     ueq_dir=zeros(m,1);
    0064     return;
    0065   end
    0066
    0067   rhs=[rhs1;rhs2];
    0068 %  dir(perm)=L' \ (D \ (L \ rhs(perm))); % this somehow gives 1x(n+m)
    0069 %  vector instead of (n+m)x1, WHY ???
    0070   dir=L' \ (D \ (L \ rhs(perm)));
    0071   dir(perm)=dir;
    0072   x_dir=dir(1:n);
    0073   ueq_dir=dir(n+1:n+m);
    0074   time_total=cputime - start_time;
    0075   nFlag=0;
    0076
    0077   obj.print(4,Inf,'LDL KKT OK, factor in %fs, total %fs, no pert=%i, pert=%.4e',time_factor, time_total, pert_try, lambda_last);
    0078   if (lambda_last>0)
    0079     obj.print(3,4,'LDL KKT solve: OK, pert=%.4e',lambda_last);
    0080   end
    0081
    0082 end
    0083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0084
    0085 function [e_pos, e_neg, e_zero] = inertia(D)
    0086 % returns the inertia triplet of matrix D, expects D to be a diagonal/block
    0087 % diagonal symmetric matrix with at most 2-times-2 blocks on the diagonal
    0088
    0089   % simple but perhaps ?slow? way how to do it
    0090   [dim dim2] = size(D);
    0091   e=eig(D);
    0092   e_pos=sum(e>0);
    0093   e_neg=sum(e<0);
    0094   e_zero=dim-e_pos-e_neg;
    0095
    0096   return;
    0097
    0098   % faster?? but more complicated way - split it to 1-by-1 / 2-by-2 blocks
    0099   e_pos=0;
    0100   e_neg=0;
    0101
    0102   i=1;
    0103   while (i<=dim)
    0104     if (i<dim && D(i,i+1))  % 2-by-2 block
    0105       e=eig(D(i:i+1,i:i+1));
    0106       e_pos = e_pos + sum(e>0);
    0107       e_neg = e_neg + sum(e<0);
    0108       i=i+2;
    0109     else           % 1-by-1 block
    0110       e_pos = e_pos + D(i,i)>0;
    0111       e_neg = e_neg + D(i,i)<0;
    0112       i=i+1;
    0113     end
    0114   end
    0115   e_zero=dim-e_pos-e_neg;
    0116
    0117 end
    0118