KerMor  0.9
Model order reduction for nonlinear dynamical systems and nonlinear approximation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MatUtils.m
Go to the documentation of this file.
1 
2 
3 /* (Autoinserted by mtoc++)
4  * This source code has been filtered by the mtoc++ executable,
5  * which generates code that can be processed by the doxygen documentation tool.
6  *
7  * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
8  * Except for the comments, the function bodies of your M-file functions are untouched.
9  * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
10  * attached source files that are highly readable by humans.
11  *
12  * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
13  * the correct locations in the source code browser.
14  * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
15  */
16 
17 class MatUtils {
32  public: /* ( Static ) */
33 
34  static function [A , idxmat ] = laplacemat(h,d1,d2) {
35 
36  m = d1*d2;
37  idxmat = zeros(d1,d2);
38  idxmat(:) = 1:m;
39 
40  inner = find(conv2(ones(size(idxmat)),[0 1 0;1 0 1; 0 1 0]," same ")==4);
41 
42  i = []; j = []; s =[];
43  /* % Inner points
44  * N E S W */
45  createStencil([-1 d1 1 -d1],[-4 1 1 1 1], inner);
46 
47  /* % Points that are on a boundaries
48  * Left boundary neumann */
49  createStencil([-1 d1 1] ,[-3 1 1 1],(2:d1-1)^t);
50  /* Right boundary neumann */
51  createStencil([-1 -d1 1] ,[-3 1 1 1],(m-d1+2:m-1)^t);
52  /* Top boundary neumann */
53  createStencil([-d1 1 d1] ,[-3 1 1 1],(d1+1:d1:m-2*d1+1)^t);
54  /* Bottom boundary neumann */
55  createStencil([-d1 -1 d1] ,[-3 1 1 1],(2*d1:d1:m-d1)^t);
56 
57  /* Edge points
58  * Top left */
59  createStencil([1 d1] ,[-2 1 1],1);
60  /* Top right */
61  createStencil([1 -d1] ,[-2 1 1],m-d1+1);
62  /* bottom left */
63  createStencil([-1 d1] ,[-2 1 1],d1);
64  /* bottom right */
65  createStencil([-1 -d1] ,[-2 1 1],m);
66  /* % Compile matrix */
67  A = (1/h^2)*sparse(i,j,s,m,m);
68 
69  function createStencil ( stencil , weights , points )
70  /*
71  % Parameters:
72  * stencil - Testdescription
73  *
74  * weights - bla bla */
75 
76  i = [i; idxmat(points)];
77  j = [j; idxmat(points)];
78  s = [s; weights(1)* ones(size(points))];
79  idx=2;
80  for offset = stencil
81  i = [i; idxmat(points)];
82  j = [j; idxmat(points+offset)];
83  s = [s; weights(idx)*ones(size(points))];
84  idx=idx+1;
85  end
86  end
87  }
104  static function matrix<double>A = divcdivumat(px,py,gradc) {
105 
106  h = diff(px);
107  if h == 0
108  error(" px must increase along first dimension. ");
109  elseif any(any(abs((h(1) - h)./h) > 1e-13)) || any(any(abs((h(1) - diff(py" ))./diff(py ")) > 1e-13))
110  error(" Need equidistant points in both directions. ");
111  end
112 
113  d1 = size(px,1);
114  d2 = size(px,2);
115  m = d1*d2;
116  idxmat = zeros(d1,d2);
117  idxmat(:) = 1:m;
118 
119  inner = find(conv2(ones(size(idxmat)),[0 1 0;1 0 1; 0 1 0]," same ")==4);
120 
121  i = []; j = []; s =[];
122  /* % Inner points
123  * N E S W */
124  createStencil([-1 d1 1 -d1], .5*[-1 1 1 -1], [1 2 1 2], inner);
125 
126  /* % Points that are on a boundaries
127  * Left boundary neumann */
128  createStencil([-1 1 0 d1] ,[-.5 .5 -1 1],[1 1 2 2],(2:d1-1)^t);
129  /* Right boundary neumann */
130  createStencil([-1 1 -d1 0] ,[-.5 .5 -1 1],[1 1 2 2],(m-d1+2:m-1)^t);
131  /* Top boundary neumann */
132  createStencil([-d1 d1 0 1] ,[-.5 .5 -1 1],[2 2 1 1],(d1+1:d1:m-2*d1+1)^t);
133  /* Bottom boundary neumann */
134  createStencil([-d1 d1 -1 0] ,[-.5 .5 -1 1],[2 2 1 1],(2*d1:d1:m-d1)^t);
135 
136  /* Edge points
137  * Top left */
138  createStencil([0 1 0 d1] ,[-1 1 -1 1],[1 1 2 2],1);
139  /* Top right */
140  createStencil([0 1 -d1 0] ,[-1 1 -1 1],[1 1 2 2],m-d1+1);
141  /* bottom left */
142  createStencil([-1 0 0 d1] ,[-1 1 -1 1],[1 1 2 2],d1);
143  /* bottom right */
144  createStencil([-1 0 -d1 0] ,[-1 1 -1 1],[1 1 2 2],m);
145  /* % Compile matrix */
146  A = sparse(i,j,s,m,m)/h(1);
147 
148  function createStencil (stencil, weights, component, points)
149  grad = zeros(length(points),2);
150  for k=1:length(points)
151  grad(k,:) = gradc([px(points(k)); py(points(k))]);
152  end
153  for pos = 1:length(stencil)
154  i = [i; idxmat(points)];/* #ok */
155 
156  j = [j; idxmat(points+stencil(pos))];/* #ok */
157 
158  s = [s; weights(pos)*abs(grad(:,component(pos)))];/* #ok */
159 
160  end
161  end
162  }
187  static function A = laplacemat3D(h,g) {
188 
189  i = []; j = []; s =[];
190  m = g.Dims(1);
191  k = g.Dims(2);
192  L = g.IndexMatrix;
193 
194  /* % Inner points */
195  createStencil([-1 m 1 -m -m*k m*k],[6 -1 -1 -1 -1 -1 -1], g.Inner);
196 
197  /* % Points that are on inner boundaries
198  * Front boundary neumann */
199  createStencil([-1 m 1 -m m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.F);
200  /* Back boundary neumann */
201  createStencil([-1 m 1 -m -m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.Ba);
202  /* Left boundary neumann */
203  createStencil([-1 1 m m*k -m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.L);
204  /* Right boundary neumann */
205  createStencil([-1 1 -m m*k -m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.R);
206  /* Top boundary neumann */
207  createStencil([1 m -m m*k -m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.T);
208  /* Bottom boundary neumann */
209  createStencil([-1 m -m m*k -m*k] ,[5 -1 -1 -1 -1 -1],g.Sides.Bo);
210 
211  /* Corner points
212  *% Front side corners
213  * Front top */
214  createStencil([1 m -m m*k] ,[4 -1 -1 -1 -1],g.Corners.FT);
215  /* Front bottom */
216  createStencil([-1 m -m m*k] ,[4 -1 -1 -1 -1],g.Corners.FBo);
217  /* Front left */
218  createStencil([1 -1 m m*k] ,[4 -1 -1 -1 -1],g.Corners.FL);
219  /* Front right */
220  createStencil([1 -1 -m m*k] ,[4 -1 -1 -1 -1],g.Corners.FR);
221 
222  /* % Back side corners
223  * Back top */
224  createStencil([1 m -m -m*k] ,[4 -1 -1 -1 -1],g.Corners.BaT);
225  /* Back bottom */
226  createStencil([-1 m -m -m*k] ,[4 -1 -1 -1 -1],g.Corners.BaBo);
227  /* Back left */
228  createStencil([1 -1 m -m*k] ,[4 -1 -1 -1 -1],g.Corners.BaL);
229  /* Back right */
230  createStencil([1 -1 -m -m*k] ,[4 -1 -1 -1 -1],g.Corners.BaR);
231 
232  /* % front-back connecting corners
233  * top left */
234  createStencil([1 m m*k -m*k] ,[4 -1 -1 -1 -1],g.Corners.TL);
235  /* top right */
236  createStencil([1 -m m*k -m*k] ,[4 -1 -1 -1 -1],g.Corners.TR);
237  /* bottom left */
238  createStencil([-1 m m*k -m*k] ,[4 -1 -1 -1 -1],g.Corners.BoL);
239  /* bottom right */
240  createStencil([-1 -m m*k -m*k] ,[4 -1 -1 -1 -1],g.Corners.BoR);
241 
242  /* % Edge points
243  * Front Top left */
244  createStencil([1 m m*k] ,[3 -1 -1 -1],g.Edges.FTL);
245  /* Front Top right */
246  createStencil([1 -m m*k] ,[3 -1 -1 -1],g.Edges.FTR);
247  /* Front bottom left */
248  createStencil([-1 m m*k] ,[3 -1 -1 -1],g.Edges.FBoL);
249  /* Front bottom right */
250  createStencil([-1 -m m*k] ,[3 -1 -1 -1],g.Edges.FBoR);
251  /* Rear Top left */
252  createStencil([1 m -m*k] ,[3 -1 -1 -1],g.Edges.BaTL);
253  /* Rear Top right */
254  createStencil([1 -m -m*k] ,[3 -1 -1 -1],g.Edges.BaTR);
255  /* Rear bottom left */
256  createStencil([-1 m -m*k] ,[3 -1 -1 -1],g.Edges.BaBoL);
257  /* Rear bottom right */
258  createStencil([-1 -m -m*k] ,[3 -1 -1 -1],g.Edges.BaBoR);
259 
260  /* % Compile matrix */
261  A = -(1/h^2)*sparse(i,j,s,g.Points,g.Points);
262 
263  /* Some comment on how createStencil works - ABOVE */
264  function createStencil ( stencil , weights , points )
265  /* Some comment on how createStencil works - BELOW
266  *
267  * Parameters:
268  * stencil - Testdescription
269  *
270  * weights - bla bla */
271 
272  i = [i; L(points)];
273  j = [j; L(points)];
274  s = [s; weights(1)* ones(size(points))];
275  cnt=2;
276  for offset = stencil
277  i = [i; L(points)]; /* #ok */
278 
279  j = [j; L(points+offset)];/* #ok */
280 
281  s = [s; weights(cnt)*ones(size(points))];/* #ok */
282 
283  cnt=cnt+1;
284  end
285  end
286  }
304  static function A = generalizedLaplacemat3D(h,g,sigma) {
305 
306  s1 = sigma(1,1)/h(1)^2; s2 = sigma(2,2)/h(2)^2; s3 = sigma(3,3)/h(3)^2;
307  /* h1 = h(1)^2; h2 = h(2)^2; h3 = h(3)^2; */
308  i = []; j = []; s =[];
309  m = g.Dims(1);
310  k = g.Dims(2);
311 
312  /* % Inner points */
313  createStencil([-1 m 1 -m -m*k m*k],[2*(s1+s2+s3) -s1 -s2 -s1 -s2 -s3 -s3], g.Inner);
314 
315  /* % Points that are on inner boundaries
316  * Front boundary neumann */
317  createStencil([-1 m 1 -m m*k] ,[2*s1+2*s2+s3 -s1 -s2 -s1 -s2 -s3],g.Sides.F);
318  /* Back boundary neumann */
319  createStencil([-1 m 1 -m -m*k] ,[2*s1+2*s2+s3 -s1 -s2 -s1 -s2 -s3],g.Sides.Ba);
320  /* Left boundary neumann */
321  createStencil([-1 1 m m*k -m*k] ,[2*s1+s2+2*s3 -s1 -s1 -s2 -s3 -s3],g.Sides.L);
322  /* Right boundary neumann */
323  createStencil([-1 1 -m m*k -m*k] ,[2*s1+s2+2*s3 -s1 -s1 -s2 -s3 -s3],g.Sides.R);
324  /* Top boundary neumann */
325  createStencil([1 m -m m*k -m*k] ,[s1+2*s2+2*s3 -s1 -s2 -s2 -s3 -s3],g.Sides.T);
326  /* Bottom boundary neumann */
327  createStencil([-1 m -m m*k -m*k] ,[s1+2*s2+2*s3 -s1 -s2 -s2 -s3 -s3],g.Sides.Bo);
328 
329  /* Corner points
330  *% Front side corners
331  * Front top */
332  createStencil([1 m -m m*k] ,[s1+2*s2+s3 -s1 -s2 -s2 -s3],g.Corners.FT);
333  /* Front bottom */
334  createStencil([-1 m -m m*k] ,[s1+2*s2+s3 -s1 -s2 -s2 -s3],g.Corners.FBo);
335  /* Front left */
336  createStencil([1 -1 m m*k] ,[2*s1+s2+s3 -s1 -s1 -s2 -s3],g.Corners.FL);
337  /* Front right */
338  createStencil([1 -1 -m m*k] ,[2*s1+s2+s3 -s1 -s1 -s2 -s3],g.Corners.FR);
339 
340  /* % Back side corners
341  * Back top */
342  createStencil([1 m -m -m*k] ,[s1+2*s2+s3 -s1 -s2 -s2 -s3],g.Corners.BaT);
343  /* Back bottom */
344  createStencil([-1 m -m -m*k] ,[s1+2*s2+s3 -s1 -s2 -s2 -s3],g.Corners.BaBo);
345  /* Back left */
346  createStencil([1 -1 m -m*k] ,[2*s1+s2+s3 -s1 -s1 -s2 -s3],g.Corners.BaL);
347  /* Back right */
348  createStencil([1 -1 -m -m*k] ,[2*s1+s2+s3 -s1 -s1 -s2 -s3],g.Corners.BaR);
349 
350  /* % front-back connecting corners
351  * top left */
352  createStencil([1 m m*k -m*k] ,[s1+s2+2*s3 -s1 -s2 -s3 -s3],g.Corners.TL);
353  /* top right */
354  createStencil([1 -m m*k -m*k] ,[s1+s2+2*s3 -s1 -s2 -s3 -s3],g.Corners.TR);
355  /* bottom left */
356  createStencil([-1 m m*k -m*k] ,[s1+s2+2*s3 -s1 -s2 -s3 -s3],g.Corners.BoL);
357  /* bottom right */
358  createStencil([-1 -m m*k -m*k] ,[s1+s2+2*s3 -s1 -s2 -s3 -s3],g.Corners.BoR);
359 
360  /* % Edge points
361  * Front Top left */
362  createStencil([1 m m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.FTL);
363  /* Front Top right */
364  createStencil([1 -m m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.FTR);
365  /* Front bottom left */
366  createStencil([-1 m m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.FBoL);
367  /* Front bottom right */
368  createStencil([-1 -m m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.FBoR);
369  /* Rear Top left */
370  createStencil([1 m -m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.BaTL);
371  /* Rear Top right */
372  createStencil([1 -m -m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.BaTR);
373  /* Rear bottom left */
374  createStencil([-1 m -m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.BaBoL);
375  /* Rear bottom right */
376  createStencil([-1 -m -m*k] ,[s1+s2+s3 -s1 -s2 -s3],g.Edges.BaBoR);
377 
378  /* % Compile matrix */
379  A = -sparse(i,j,s,g.Points,g.Points);
380 
381  /* Some comment on how createStencil works - ABOVE */
382  function createStencil ( stencil , weights , points )
383  /* Some comment on how createStencil works - BELOW
384  *
385  * Parameters:
386  * stencil - Testdescription
387  *
388  * weights - bla bla */
389 
390  i = [i; points];
391  j = [j; points];
392  s = [s; weights(1)* ones(size(points))];
393  cnt=2;
394  for offset = stencil
395  i = [i; points]; /* #ok */
396 
397  j = [j; points+offset];/* #ok */
398 
399  s = [s; weights(cnt)*ones(size(points))];/* #ok */
400 
401  cnt=cnt+1;
402  end
403  end
404  }
424  static function Kinv = getExtendedInverse(Kinv,v,mode) {
425  if nargin < 3
426  /* Experiments show that mode 3 is best for gaussian kernel matrices */
427  mode = 3;
428  end
429  oldn = size(v,1)-1;
430  vold = v(1:oldn);
431  e = zeros(size(v));
432  e(end) = 1;
433  v(end) = (v(end)-1)/2;
434  U = [v e];
435  V = [e v]^t;
436  tmp = [Kinv zeros(oldn,1); zeros(1,oldn) 1];
437 
438  kh = vold^t*Kinv*vold;
439  det = (1 + v(end) - kh);
440  if mode == 1
441  hlp = [-vold*vold" vold; vold " 2*v(end)-vold^t*Kinv*vold]/det;
442  Kinv = tmp*(eye(oldn+1)-hlp*tmp);
443  return;
444  end
445 
446  b = (1/det) * [1+v(end) -1; -kh-v(end)^2 1+v(end)];
447  if mode == 2
448  Kinv = tmp - (tmp*U)*b*(V*tmp);
449  return;
450  end
451 
452  if mode == 3
453  a = [Kinv*vold zeros(oldn,1); v(end) 1];
454  c = [zeros(1,oldn) 1; vold^t*Kinv v(end)];
455  Kinv = tmp - a*b*c;
456  end
457  }
467  static function Asparse = toSparse(matrix<double> A,rowvec<integer> rowidx,n) {
468  o = size(A,2);
469  i = repmat(rowidx,o,1);
470  j = reshape(repmat(1:o,length(rowidx),1),[],1);
471  Asparse = sparse(i,j,A(:),n,o);
472  }
485  public: /* ( Static ) */
486 
487 
488  static function res = test_DivCDivUMat() {
489  [x,y] = meshgrid(-1:1);
490  x = x" ; y = y ";
491  gradc = @(x)[1 1];
492 
493  A = MatUtils.divcdivumat(x,y,gradc);
494 
495  gradc = @(x)[0 1];
496 
497  A1 = MatUtils.divcdivumat(x,y,gradc);
498 
499  gradc = @(x)[1 0];
500 
501  A2 = MatUtils.divcdivumat(x,y,gradc);
502 
503  res = isequal(A, A1 + A2);
504  }
505 
506 
507  static function test_ExtendedInverseDirect() {
508 
509  /* Matrix size */
510  s = 100;
511 
512  /* Experiments */
513  n = 30;
514 
515  /* Kernel widths */
516  w = 20;
517 
518  g = kernels.GaussKernel();
519  gammas = linspace(sqrt(s)/10,3*sqrt(s),w);
520  err = zeros(3,n*w);
521  for gam=1:w
522  g.Gamma= gammas(gam);
523  /* fprintf('Using gamma = %f...\n',g.Gamma); */
524  for exp=1:n
525  X = rand(s,s);
526  A = g.evaluate(X, X);
527  y = rand(s,1);
528  v = g.evaluate(X, y);
529  v(end+1) = g.evaluate(y,y);/* #ok */
530 
531  Bi = inv([A v(1:end-1); v^t]);
532  i = (gam-1)*n + exp;
533  Ai = inv(A);
534  err(1,i) = errfun(Bi,MatUtils.getExtendedInverse(Ai, v, 1));
535  err(2,i) = errfun(Bi,MatUtils.getExtendedInverse(Ai, v, 2));
536  err(3,i) = errfun(Bi,MatUtils.getExtendedInverse(Ai, v, 3));
537  end
538  end
539  e1best = sum(err(1,:) < min(err(2:3,:),[],1));
540  e2best = sum(err(2,:) < min(err([1 3],:),[],1));
541  e3best = sum(err(3,:) < min(err(1:2,:),[],1));
542  semilogy(1:n*w,err);
543  title(sprintf(" Plot for %d gamma values with matrix size %d, %d experiments each.\nBest inverse (norm-difference wise): 1:%d, 2:%d, 3:%d ",w,s,n,e1best,e2best,e3best));
544  axis tight;
545 
546  function e = errfun(A,B)
547  /* e = norm(A-B); */
548  e = max(abs(A(:)-B(:)));
549  end
550  }
560  static function test_ExtendedInverseSequential() {
561 
562  /* Matrix size */
563  s = 200;
564 
565  /* Kernel widths */
566  w = 6;
567 
568  g = kernels.GaussKernel();
569  gammas = linspace(sqrt(s)/10,4*sqrt(s),w);
570  err = zeros(3,w*(s-1));
571  T = zeros(4,w*(s-1));
572  X = rand(s,s);
573  for gam=1:w
574  g.Gamma= gammas(gam);
575  /* fprintf('Using gamma = %f...\n',g.Gamma); */
576  Am = g.evaluate(X,X);
577 /* fx = rand(s,1); */
578  Ai1 = 1/Am(1,1);
579  Ai2 = 1/Am(1,1);
580  Ai3 = 1/Am(1,1);
581  for part=1:s-1
582  i = (gam-1)*(s-1) + part;
583  /* A = Am(1:part,1:part); */
584  v = Am(1:part+1,part+1);
585  t = tic;
586  Ai1 = MatUtils.getExtendedInverse(Ai1, v, 1);
587  T(1,i) = toc(t);
588  t = tic;
589  Ai2 = MatUtils.getExtendedInverse(Ai2, v, 2);
590  T(2,i) = toc(t);
591  t = tic;
592  Ai3 = MatUtils.getExtendedInverse(Ai3, v, 3);
593  T(3,i) = toc(t);
594  t = tic;
595  Bi = inv(Am(1:part+1,1:part+1));
596  T(4,i) = toc(t);
597 
598  err(1,i) = errfun(Bi,Ai1);
599  err(2,i) = errfun(Bi,Ai2);
600  err(3,i) = errfun(Bi,Ai3);
601  end
602  end
603  e1best = sum(err(1,:) < min(err(2:3,:),[],1));
604  e2best = sum(err(2,:) < min(err([1 3],:),[],1));
605  e3best = sum(err(3,:) < min(err(1:2,:),[],1));
606  semilogy(1:w*(s-1),err);
607  T = sum(T,2);
608  title(sprintf([" Plot for %d gamma values with matrix size %d.\n "...
609  " Count for best mode (norm-difference wise): 1:%d, 2:%d, 3:%d\n "...
610  " Times for modes: 1: %1.3fs, 2: %1.3fs, 3: %1.3fs, direct: %1.3fs "],w,s,e1best,e2best,e3best,...
611  T(1),T(2),T(3),T(4)));
612  axis tight;
613 
614 
615  function e = errfun(A,B)
616  /* e = norm(A-B); */
617  e = max(abs(A(:)-B(:)));
618  end
619  }
620 
621 
622 
623 };
624 
625 
626 
MatUtils.toSparse
static function Asparse = toSparse(matrix< double > A,rowvec< integer > rowidx, n)
Converts a full matrix A to a sparse matrix, where the entries of A are split up to the row indices s...
Definition: MatUtils.m:467
MatUtils.test_ExtendedInverseDirect
static function test_ExtendedInverseDirect()
Tests direct inversion extension for a matrix of size s for a fixed number of experiments each...
Definition: MatUtils.m:507
MatUtils.test_ExtendedInverseSequential
static function test_ExtendedInverseSequential()
Definition: MatUtils.m:560
MatUtils.test_DivCDivUMat
static function res = test_DivCDivUMat()
Definition: MatUtils.m:488
MatUtils.laplacemat
static function [ A , idxmat ] = laplacemat(h, d1, d2)
Computes a 2D diffusion sparse matrix with zero neuman boundary conditions.
Definition: MatUtils.m:34
X
#define X(i, j)
Definition: dontuse_evaluate.c:5
MatUtils.generalizedLaplacemat3D
static function A = generalizedLaplacemat3D(h, g, sigma)
Computes a 3D generalized laplacian sparse matrix .
Definition: MatUtils.m:304
MatUtils
MatUtils: Matrix utility functions.
Definition: MatUtils.m:17
MatUtils.getExtendedInverse
static function Kinv = getExtendedInverse(Kinv, v, mode)
For a matrix K with known inverse Kinv the matrix is extended by the vector.
Definition: MatUtils.m:424
k
*for k
Definition: Gordon66SarcoForceLength.m:114
t
* t
Definition: Gordon66SarcoForceLength.m:116
MatUtils.divcdivumat
static function matrix< double > A = divcdivumat(px, py, gradc)
Computes the 2D matrix for the expression arising in the general spatial-dependent diffusion term ...
Definition: MatUtils.m:104
MatUtils.laplacemat3D
static function A = laplacemat3D(h, g)
Computes a 3D laplacian sparse matrix.
Definition: MatUtils.m:187