not computed L yet

src/LAPACK/LUdecomposition.winxed

@@ -0,0 +1,135 @@
+$load "inverse.pbc";
+
+namespace ludecompo_func{
+
+const int PRINT_DEBUG_STUFF = 0;
+
+function Uform(var A,var u)
+{
+	var pla = loadlib("linalg_group");
+	int m,n;
+	m=A.rows;
+	n=A.cols;
+	//copy(A,u,m,n);	
+	u=A.clone();
+	int i,j,k,pivoitrow;
+	var temp;
+	var index_temp;
+	for(i=0;i<m;++i)
+	{
+		pivoitrow=i;
+		for(j=i+1;j<m;++j)
+		{
+			if(abs(u[i,j])>abs(u[pivoitrow,i]))
+				pivoitrow=j;
+		}
+		for(k=i;k<n;++k)
+		{
+			temp=u[i,k];
+			u[i,k]=u[pivoitrow,k];
+			u[pivoitrow,k]=temp;
+		}
+		for(j=i+1;j<m;++j)
+		{
+			index_temp=u[j,i]/u[i,i];
+			for(k=i;k<n;++k)
+			{
+				u[j,k]=u[j,k]-(u[i,k]*index_temp);
+			}
+		}
+	}
+
+}
+
+
+function Lform(var A,var l,var u)
+{
+        var pla = loadlib("linalg_group");
+        int m,n;
+        m=A.rows;
+        n=A.cols;
+	copy(A,l,m,n);
+
+        int i,j,k,pivoitrow;
+        var temp;
+        var index_temp;
+/*        for(i=m-1;i>0;--i)
+        {
+                pivoitrow=i;
+		say(pivoitrow);
+                for(j=n-1-(m-i-1);j>=0;--j)
+                {
+                        if(abs(l[i,j])>abs(l[pivoitrow,i]))
+                                pivoitrow=j;
+                }
+		say(pivoitrow);
+                for(k=0;k<n;++k)
+                {
+                        temp=l[i,k];
+                        l[i,k]=l[pivoitrow,k];
+                        l[pivoitrow,k]=temp;
+                }
+                for(j=i-1;j>=0;--j)
+                {
+                        index_temp=l[j,i]/l[i,i];
+                        for(k=0;k<i;++k)
+                        {
+                                l[j,k]=l[j,k]-(l[i,k]*index_temp);
+                        }
+                }
+        }
+*/
+
+	using dgetri_func.dgetri_exec;
+        int info;
+        info=dgetri_exec(u);
+	l=A*u;
+
+
+}
+
+function abs(var a)
+{
+	if(a<0)
+		return (-1*a);
+	return (a);
+}
+
+function copy(var a,var u,var m,var n)
+{
+	int i,j;
+	for(i=0;i<m;++i)
+		for(j=0;j<n;++j)
+			u[i,j]=a[i,j];
+}
+
+}
+
+
+
+
+function main[main](var args)
+{
+
+        var pla = loadlib("linalg_group");
+        var a = new 'NumMatrix2D';
+        var u=new 'NumMatrix2D';
+        var l=new 'NumMatrix2D';
+	a.initialize_from_args(3, 3,
+                                2.0, -1.0, 1.0,
+                                4.0, 1.0, -1.0,
+                                1.0, 1.0, 1.0);
+	using ludecompo_func.Uform;
+	Uform(a,u);
+	say("Successful");
+	say("U=");
+	say(u);
+	//say("A=");
+	//say(a);
+	using ludecompo_func.Lform;
+        Lform(a,l,u);
+        say("Successful");
+        say("L=");
+        say(l);
+
+}