1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
## Copyright (C) 2006,2007,2008 Carlo de Falco, Culpo Massimiliano
##
## This file is part of:
## OCS - A Circuit Simulator for Octave
##
## OCS is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program (see the file LICENSE); if not,
## see <http://www.gnu.org/licenses/>.
##
## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
## author: culpo@math.uni-wuppertal.de
## -*- texinfo -*-
## @deftypefn{Function File} {[@var{A},@var{Jac},@var{res}] =}@
## asm_build_system(@var{instruct},@var{x},@var{t})
##
## Cycle through the circuit description structure @var{instruct}
## to build the system matrices @var{A}, @var{Jac}, @var{res} for
## the current step of the Newton method.
##
## @var{x} is the current value of the state variables, while @var{t} is the current time point.
##
## See the @cite{IFF file format specifications} for details about
## the output matrices.
##
## @seealso{asm_initialize_system,prs_iff}
##
## @end deftypefn
function [A,Jac,res] = asm_build_system(instruct,x,t);
## Check input
if nargin != 3
error("asm_build_system: wrong number of input parameters.");
elseif !(isstruct(instruct) && isfield(instruct,"LCR") &&
isfield(instruct,"NLC"))
error("asm_build_system: first input is not a valid structure.");
elseif !isvector(x)
error("asm_build_system: second input is not a valid vector.");
elseif !isscalar(t)
error("asm_build_system: third input is not a valid scalar.");
endif
n = instruct.totextvar + instruct.totintvar;
A = sparse(n,n);
Jac = sparse(n,n);
res = sparse(n,1);
## NLC section
nblocks = length(instruct.NLC);
for ibl = 1:nblocks
for iel = 1:instruct.NLC(ibl).nrows
## Evaluate element
if instruct.NLC(ibl).nintvar(iel)
intvars = instruct.totextvar+instruct.NLC(ibl).osintvar(iel) + ...
[1:instruct.NLC(ibl).nintvar(iel)]';
else
intvars=[];
endif
il = instruct.NLC(ibl).vnmatrix(iel,:)';
nzil = find(il!=0);
y = zeros(size(il));
y(nzil) = x(il(nzil));
z = x(intvars);
[a,b,c] = feval(instruct.NLC(ibl).func,...
instruct.NLC(ibl).section,...
instruct.NLC(ibl).pvmatrix(iel,:),...
instruct.NLC(ibl).parnames,...
y,z,t);
## Assemble matrices
## Global indexing
vars = [il(nzil);intvars];
## Local indexing
lclvars = [nzil; instruct.NLC(ibl).nextvar + (1:length(intvars))' ];
## Reshaping sparse stamps
a = a(lclvars,lclvars);
b = b(lclvars,lclvars);
c = reshape(c(lclvars),[],1);
[na,ma,va] = find(a);
[nb,mb,vb] = find(b);
[nc,mc,vc] = find(c);
## Stamping
A += sparse(vars(na),vars(ma),va,n,n);
Jac += sparse(vars(nb),vars(mb),vb,n,n);
res += sparse(vars(nc),1,vc,n,1);
endfor
endfor
endfunction
|
