[Octave-cvsupdate] octave-forge/main/splines Makefile,1.4,1.5 csape.m,1.4,1.5 dgtsv.f,1.1.1.1,NONE dptsv.f,1.1.1.1,NONE dpttrf.f,1.1.1.1,NONE dpttrs.f,1.1.1.1,NONE dptts2.f,1.1.1.1,NONE trisolve.cc,1.6,NONE trisolve.m,1.2,NONE trisolve.tst,1.1.1.1,NONE

[David B. <ad...@us...>]

Update of /cvsroot/octave/octave-forge/main/splines
In directory sc8-pr-cvs3.sourceforge.net:/tmp/cvs-serv9045

Modified Files:
	Makefile csape.m 
Removed Files:
	dgtsv.f dptsv.f dpttrf.f dpttrs.f dptts2.f trisolve.cc 
	trisolve.m trisolve.tst 
Log Message:
Remove trisolve and dependency on trisolve and use sparse tridiagonal solver in octave 2.9.x instead

--- dpttrs.f DELETED ---

Index: csape.m
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/splines/csape.m,v
retrieving revision 1.4
retrieving revision 1.5
diff -u -d -r1.4 -r1.5
--- csape.m	5 Feb 2004 05:17:57 -0000	1.4
+++ csape.m	25 May 2006 19:50:44 -0000	1.5
@@ -80,7 +80,8 @@
       g(n - 2,:) = 3 / 2 * (3 * (a(n,:) - a(n - 1,:)) / h(n - 1) - valc(2))\
           - 3 * (a(n - 1,:) - a(n - 2,:)) / h(n - 2);
 
-      c(2:n - 1,:) = trisolve(dg,e,g);
+      c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
+
     end
 
     c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * valc(1) 
@@ -114,7 +115,7 @@
     else
       dg = 2 * (h(1:n - 2) .+ h(2:n - 1));
       e = h(2:n - 2);
-      c(2:n - 1,:) = trisolve (dg,e,g);
+      c(2:n - 1,:) = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-2,n-2) \ g;
     end
         
     b(1:n - 1,:) = diff (a) ./ h(1:n - 1,idx)\
@@ -136,7 +137,18 @@
     if (n > 3)
       dg = 2 * (h(1:n - 1) .+ h(2:n));
       e = h(2:n - 1);
-      c(2:n,idx) = trisolve(dg,e,g,h(1),h(1));
+
+      ## Use Sherman-Morrison formula to extend the solution
+      ## to the cyclic system. See Numerical Recipes in C, pp 73-75
+      gamma = - dg(1);
+      dg(1) -=  gamma;
+      dg(end) -= h(1) * h(1) / gamma; 
+      z = spdiags([[e(:);0],dg,[0;e(:)]],[-1,0,1],n-1,n-1) \ ...
+	  [[gamma; zeros(n-3,1); h(1)],g];
+      fact = (z(1,2:end) + h(1) * z(end,2:end) / gamma) / ...
+	  (1.0 + z(1,1) + h(1) * z(end,1) / gamma);
+
+      c(2:n,idx) = z(:,2:end) - z(:,1) * fact;
     endif
 
     c(1,:) = c(n,:);
@@ -167,14 +179,14 @@
       ldg = udg = h(2:n - 2);
       udg(1) = udg(1) - h(1);
       ldg(n - 3) = ldg(n-3) - h(n - 1);
-      c(2:n - 1,:) = trisolve(ldg,dg,udg,g);
+      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
 
     elseif (n == 4)
 
       dg = [h(1) + 2 * h(2), 2 * h(2) + h(3)];
       ldg = h(2) - h(3);
       udg = h(2) - h(1);
-      c(2:n - 1,:) = trisolve(ldg,dg,udg,g);
+      c(2:n - 1,:) = spdiags([[ldg(:);0],dg,[0;udg(:)]],[-1,0,1],n-2,n-2) \ g;
       
     else # n == 3
 	    

--- trisolve.tst DELETED ---

--- dptsv.f DELETED ---

--- dgtsv.f DELETED ---

--- trisolve.cc DELETED ---

Index: Makefile
===================================================================
RCS file: /cvsroot/octave/octave-forge/main/splines/Makefile,v
retrieving revision 1.4
retrieving revision 1.5
diff -u -d -r1.4 -r1.5
--- Makefile	7 Jul 2004 09:38:49 -0000	1.4
+++ Makefile	25 May 2006 19:50:44 -0000	1.5
@@ -1,12 +1,8 @@
 sinclude ../../Makeconf
 
-TRISOLVE_OBJ = trisolve.o dgtsv.o dptsv.o dpttrf.o dpttrs.o dptts2.o
 PCHIP_OBJ = pchip_deriv.o dpchim.o dpchst.o
 
-all: trisolve.oct pchip_deriv.oct
-
-trisolve.oct: $(TRISOLVE_OBJ)
-	$(MKOCTFILE) -v -o $@ $(TRISOLVE_OBJ)
+all: pchip_deriv.oct
 
 pchip_deriv.oct: $(PCHIP_OBJ)
 	$(MKOCTFILE) -v -o $@ $(PCHIP_OBJ)

--- dpttrf.f DELETED ---

--- trisolve.m DELETED ---

--- dptts2.f DELETED ---