← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/algebre/chfact.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/algebre/chfact.sci
 
1
function spcho=chfact(A)
 
2
//cholesky factors, returned in a tlist
 
3
//spcho  = {xlnz, nnzl, xsuper, xlindx, lindx, snode,
 
4
//          split, tmpsiz, perm, invp, lnz}.
 
5
//
 
6
//  invp=spcho('invp');xlnz=spcho('xlnz');
 
7
//  xlindx=spcho('xlindx');lindx=spcho('lindx');lnz=spcho('lnz');
 
8
//  P=sparse([(1:m)',invp],ones(invp),[m,m]);
 
9
//  adjncy = spcompack(xlnz,xlindx,lindx);
 
10
//  R=adj2sp(xlnz,adjncy,lnz);
 
11
//  =>P*R*R'*P'=A
 
12
m = size(A,1);
 
13
if size(A,2) ~= m | type(A)~=5 | A == [] then,
 
14
   error('Matrix must be square, sparse and nonempty.');
 
15
 end;
 
16
neqns=m;
 
17
[xadj,adjncy,v]=sp2adj(A-diag(diag(A)));
 
18
[perm,invp,nofsub]=ordmmd(xadj,adjncy,neqns);
 
19
cachsz = 16;
 
20
spcho=symfct(xadj,adjncy,perm,invp,cachsz,neqns);
 
21
[xadjf,adjncyf,v]=sp2adj(A);
 
22
spcho=inpnv(xadjf,adjncyf,v,spcho);
 
23
level=8;
 
24
spcho=blkfc1(spcho,level);
 
25
 
 
26
function [spcho]= blkfc1(spcho,level)
 
27
//retrieves Fortran variables (see sfinit.f,bfinit.f,symfct.f )
 
28
//[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12);
 
29
xsuper=spcho('xsuper');
 
30
nsuper=size(xsuper,1)-1;
 
31
snode=spcho('snode');
 
32
neqns=size(snode,1);
 
33
lnz=spcho('lnz');
 
34
nnzl=size(lnz,1);
 
35
iwsiz = 2*neqns+2*nsuper;
 
36
iwork = zeros(iwsiz,1);
 
37
tmpsiz = spcho('tmpsiz');
 
38
tmpvec= zeros(tmpsiz,1);
 
39
iflag=0;
 
40
// calling blkfc1i primitive
 
41
[lnz,iflag]=blkfc1i(neqns,nsuper,xsuper,snode,spcho('split'),spcho('xlindx'),spcho('lindx'),spcho('xlnz'),lnz,iwsiz,iwork,tmpsiz,tmpvec,iflag,level);
 
42
//
 
43
if max(abs(lnz)) > 5d63 then
 
44
   warning('  Possible matrix is not positive definite');
 
45
end;
 
46
 
 
47
select iflag
 
48
case -1 then, 
 
49
  error('nonpositive diag. encountered, matrix is not positive def'),
 
50
case -2 then, 
 
51
  error('Insufficient working storage in blkfct, temp(*)'),
 
52
case -3 then, 
 
53
  error('Insufficient working storage in blkfct, iwork(*)'),
 
54
end;
 
55
//
 
56
spcho('lnz')=lnz;
 
57
 
 
58
function rhs=blkslv(spcho,rhs)
 
59
//
 
60
//[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12);
 
61
xsuper=spcho('xsuper');
 
62
nsuper=size(xsuper,1)-1;
 
63
neqns =size(rhs,1);
 
64
rhs=blkslvi(nsuper,xsuper,spcho('xlindx'),spcho('lindx'),spcho('xlnz'),...
 
65
spcho('lnz'),rhs);
 
66
 
 
67
 
 
68
function [spcho]=inpnv(xadjf,adjf,anzf,spcho)
 
69
//
 
70
//[xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,lnz]=spcho(2:12);
 
71
//
 
72
xsuper=spcho('xsuper');
 
73
neqns=size(xadjf,1)-1,
 
74
nsuper=size(xsuper,1)-1,
 
75
//
 
76
offset=zeros(neqns,1);
 
77
lnz=zeros(spcho('nnzl'),1); 
 
78
//
 
79
lnz = inpnvi(neqns,xadjf,adjf,anzf,spcho('perm'),spcho('invp'),nsuper,...
 
80
xsuper,spcho('xlindx'),spcho('lindx'),spcho('xlnz'),lnz,offset);
 
81
//
 
82
spcho('lnz')=lnz;
 
83
 
 
84
function [spcho] = symfct(xadj,adjncy,perm,invp,cachsz,neqns)
 
85
//
 
86
// sfinit - input
 
87
//
 
88
nnza=size(adjncy,1);
 
89
iwsiz  = 7*neqns+4;
 
90
iwork=zeros(iwsiz,1);
 
91
///
 
92
if size(perm)~= [neqns,1] then, error(' SYMFCT requires PERM to be neqns x 1'),
 
93
end;
 
94
if size(invp)~= [neqns,1] then, error(' SYMFCT requires INVN to be neqns x 1'),
 
95
end;
 
96
if size(cachsz)~= [1,1] then, error(' SYMFCT requires CACHSZ  to be 1 x 1'),
 
97
end;
 
98
//
 
99
[perm,invp,colcnt,nnzl,nsub,nsuper,snode,xsuper,iflag]=...
 
100
sfinit(neqns,nnza,xadj,adjncy,perm,invp,iwsiz,iwork);
 
101
//
 
102
if iflag == -1 then error(' Insufficient working storage in sfinit'),end;
 
103
//
 
104
bb=xsuper(1:nsuper+1,1);
 
105
xsuper=bb
 
106
iwsiz  = 2*nsuper+2*neqns+1;
 
107
iwork=zeros(iwsiz,1);
 
108
//
 
109
//
 
110
[xlindx,lindx,xlnz,iflag]=symfcti(neqns,nnza,xadj,adjncy,perm,...
 
111
invp,colcnt,nsuper,xsuper,snode,nsub,iwsiz,iwork)
 
112
if iflag == -1 then error(' Insufficient working storage in symfct'),end;
 
113
if iflag == -2 then error(' Inconsistancy in the input in symfct'),end;
 
114
//
 
115
[tmpsiz,split]=bfinit(neqns,nsuper,xsuper,snode,xlindx,lindx,cachsz)
 
116
 
 
117
spcho = tlist(['chol','xlnz','nnzl','xsuper','xlindx','lindx','snode','split',...
 
118
'tmpsiz','perm','invp','lnz'],... 
 
119
xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp,[])