octave_packages/ocs-0.1.3/asm/asm_build_system.m

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