← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to demos/optimization/lmitool/bigM.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
demos/optimization/lmitool/bigM.sci
 
1
function [x,Z,z,ul,iters]=bigM(F,blck_szs,c,x0,M,nu,abstol,reltol,tv,maxiters);
 
2
// [x,Z,z,ul,iters]=bigM(F,blck_szs,c,x0,M,nu,abstol,reltol,tv,maxiters);
 
3
//
 
4
// minimize    c^T x 
 
5
// subject to  F(x) = F0 + x1*F1 + ... + xm*Fm >= 0 
 
6
//             Tr F(x) <= M
 
7
//
 
8
// maximize    -Tr F0*(Z-zI) - Mz
 
9
// subject to  Tr Fi*(Z-zI) = c_i
 
10
//             Z >= 0, z>= 0
 
11
//
 
12
// Convergence criteria:
 
13
// (1) maxiters is exceeded
 
14
// (2) duality gap is less than abstol
 
15
// (3) primal and dual objective are both positive and
 
16
//     duality gap is less than (reltol * dual objective)
 
17
//     or primal and dual objective are both negative and
 
18
//     duality gap is less than (reltol * minus the primal objective)
 
19
// (4) reltol is negative and
 
20
//     primal objective is less than tv or dual objective is greater
 
21
//     than tv
 
22
//
 
23
// Input arguments:
 
24
// F:         (sum_i n_i^2) times (m+1) matrix
 
25
//            [ F_0^1(:) F_1^1(:) ... F_m^1(:) ]
 
26
//            [ F_0^2(:) F_1^2(:) ... F_m^2(:) ]
 
27
//                ...      ...          ...
 
28
//            [ F_0^L(:) F_1^L(:) ... F_m^L(:) ]
 
29
//            F_i^j: jth block of F_i, size n_i times n_i.
 
30
// blck_szs:  L-vector [n_1 ... n_L], dimensions of diagonal blocks.
 
31
// c:         m-vector.  Specifies primal objective.
 
32
// x0:        m-vector.  The primal starting point.  F(x0) > 0.  
 
33
// M:         scalar. M > Tr F(x0).    
 
34
// nu:        >= 1.0.  Controls the rate of convergence.
 
35
// abstol:    absolute tolerance.
 
36
// reltol:    relative tolerance.  Has a special meaning when negative.
 
37
// tv:        target value.
 
38
// maxiters:  maximum number of iterations.
 
39
//
 
40
// Output arguments:
 
41
// x:         m-vector; last primal iterate.
 
42
// Z:         last dual iterate; block-diagonal matrix stored as 
 
43
//            [ Z^1(:);  Z^2(:); ... ; Z^L(:) ].
 
44
// z:         scalar part of last dual iterate.  
 
45
// ul:        ul(1): primal objective, ul(1): dual objective.
 
46
// iters:     number of iterations taken.
 
47
 
 
48
// Copyright INRIA
 
49
 
 
50
[rowf,colf]=size(F);
 
51
m = colf-1;
 
52
if (rowf ~= sum(blck_szs.*blck_szs))
 
53
    error('Dimensions of F do not match blck_szs.');
 
54
end;
 
55
 
 
56
[rowx0,colx0]=size(x0);
 
57
if (rowx0 ~= m) | (colx0 ~= 1)
 
58
   error('x0 must be an m-vector.'); 
 
59
end;
 
60
 
 
61
if (prod(size(x0)) ~= m), 
 
62
   error('c must be an m-vector.'); 
 
63
end;
 
64
 
 
65
// I is the identity
 
66
I = zeros(rowf,1);
 
67
blck_szs=matrix(blck_szs,1,prod(size(blck_szs)));
 
68
k=0;  for n=blck_szs,
 
69
   I(k+[1:n*n]) = matrix(eye(n,n),n*n,1);   // identity
 
70
   k = k+n*n;   // k = sum n_i*n_i
 
71
end;
 
72
 
 
73
// Z0 = projection of I on dual feasible space 
 
74
 
 
75
Z0 = I-F(:,2:m+1) * ...
 
76
     ( (F(:,2:m+1)'*F(:,2:m+1)) \ ( F(:,2:m+1)'*I - c ) );
 
77
 
 
78
// mineigZ is the smallest eigenvalue of Z0
 
79
mineigZ = 0.0;
 
80
k=0; for n=blck_szs,
 
81
  mineigZ = min(mineigZ, min(real(spec(matrix(Z0(k+[1:n*n]),n,n)))));
 
82
  k=k+n*n;
 
83
end;
 
84
 
 
85
// z = max( 1e-5, -1.1*mineigZ )
 
86
Z0(k+1) = max( 1e-5, -1.1*mineigZ);  
 
87
Z0(1:k) = Z0(1:k) + Z0(k+1)*I; 
 
88
 
 
89
if (M < I'*F*[1;x0] + 1e-5), 
 
90
   error('M must be strictly greater than trace of F(x0).'); 
 
91
end;
 
92
 
 
93
// add scalar block Tr F(x) <= M
 
94
F = [F; M-I'*F(:,1),-I'*F(:,2:m+1)]; 
 
95
blck_szs = [blck_szs,1];
 
96
[x,Z,ul,info]=...
 
97
  semidef(x0,pack(Z0),pack(F),blck_szs,c,[nu,abstol,reltol,tv,maxiters]);
 
98
iters = info(2);
 
99
nz=prod(size(Z))
 
100
z=Z(nz)
 
101
Z=unpack(Z(1:nz-1),blck_szs(1:prod(size(blck_szs))-1))
 
102
Z = Z(1:k);
 
103
 
 
104
 
 
105