[Scipy-svn] r3253 - in trunk/scipy/sandbox/multigrid: . multigridtools tests

scipy-svn@scip... scipy-svn@scip...
Wed Aug 22 10:55:17 CDT 2007


Author: wnbell
Date: 2007-08-22 10:55:13 -0500 (Wed, 22 Aug 2007)
New Revision: 3253

Modified:
   trunk/scipy/sandbox/multigrid/
   trunk/scipy/sandbox/multigrid/multigrid.py
   trunk/scipy/sandbox/multigrid/multigridtools/
   trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i
   trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/simple_test.py
   trunk/scipy/sandbox/multigrid/tests/
Log:
changed properties




Property changes on: trunk/scipy/sandbox/multigrid
___________________________________________________________________
Name: svn:ignore
   + *.so
*.bak
*.pyc
ignore.txt



Modified: trunk/scipy/sandbox/multigrid/multigrid.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multigrid.py	2007-08-21 22:12:51 UTC (rev 3252)
+++ trunk/scipy/sandbox/multigrid/multigrid.py	2007-08-22 15:55:13 UTC (rev 3253)
@@ -8,6 +8,19 @@
     
 from pydec import gauss_seidel,diag_sparse,inf_norm
 
+
+def poisson_problem1D(N):
+    """
+    Return a sparse CSC matrix for the 2d N*N poisson problem
+    with standard 5-point finite difference stencil
+    """
+    D = 2*numpy.ones(N)
+    O =  -numpy.ones(N)
+    return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N)
+
+
+
+
 def poisson_problem(N):
     """
     Return a sparse CSC matrix for the 2d N*N poisson problem
@@ -51,6 +64,17 @@
     return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
 
 
+def sa_no_threshold(A):
+    if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
+    
+    #tentative (non-smooth) interpolation operator I
+    Ij = multigridtools.sa_get_aggregates(A.shape[0],A.indptr,A.indices)
+    Ip = numpy.arange(len(Ij)+1)
+    Ix = numpy.ones(len(Ij))
+    
+    return scipy.sparse.csr_matrix((Ix,Ij,Ip))
+
+
 def sa_constant_interpolation(A,epsilon=0.08):
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
     


Property changes on: trunk/scipy/sandbox/multigrid/multigridtools
___________________________________________________________________
Name: svn:ignore
   + *.so
*.bak
*.pyc
ignore.txt



Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2007-08-21 22:12:51 UTC (rev 3252)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools.i	2007-08-22 15:55:13 UTC (rev 3253)
@@ -16,7 +16,7 @@
 
 %feature("autodoc", "1");
 
-%include "../../../sparse/sparsetools/numpy.i"
+%include "numpy.i"
 
 %init %{
     import_array();
@@ -61,7 +61,7 @@
 %enddef
 
 
-I_IN_ARRAY1( int  )
+I_IN_ARRAY1( int         )
 T_IN_ARRAY1( float       )
 T_IN_ARRAY1( double      )
 
@@ -70,8 +70,7 @@
  /*
   * OUT types
   */
-%define I_ARRAY_ARGOUT( ctype, atype )
-VEC_ARRAY_ARGOUT( ctype, atype )
+%define I_ARRAY_ARGOUT( ctype )
 %apply std::vector<ctype>* array_argout {
     std::vector<ctype>* Ap,
     std::vector<ctype>* Ai,
@@ -91,8 +90,7 @@
 };
 %enddef
 
-%define T_ARRAY_ARGOUT( ctype, atype )
-VEC_ARRAY_ARGOUT( ctype, atype )
+%define T_ARRAY_ARGOUT( ctype )
 %apply std::vector<ctype>* array_argout {
     std::vector<ctype>* Ax, 
     std::vector<ctype>* Bx,
@@ -104,9 +102,9 @@
 };
 %enddef
 
-I_ARRAY_ARGOUT( int,   INT)
-T_ARRAY_ARGOUT( float,       FLOAT   )
-T_ARRAY_ARGOUT( double,      DOUBLE  )
+I_ARRAY_ARGOUT( int    )
+T_ARRAY_ARGOUT( float  )
+T_ARRAY_ARGOUT( double )
 
 
 

Modified: trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2007-08-21 22:12:51 UTC (rev 3252)
+++ trunk/scipy/sandbox/multigrid/multigridtools/multigridtools_wrap.cxx	2007-08-22 15:55:13 UTC (rev 3253)
@@ -2576,7 +2576,9 @@
 #endif
 #include "stdio.h"
 #include <numpy/arrayobject.h>
+#include "complex_ops.h"
 
+
 /* The following code originally appeared in enthought/kiva/agg/src/numeric.i,
  * author unknown.  It was translated from C++ to C by John Hunter.  Bill
  * Spotz has modified it slightly to fix some minor bugs, add some comments
@@ -3089,7 +3091,7 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg4)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg4))[0]),sizeof(int)*length);	 
     delete arg4; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -3192,21 +3194,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg6))[0]),sizeof(int)*length);	 
     delete arg6; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg7)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg7))[0]),sizeof(int)*length);	 
     delete arg7; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg8)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_FLOAT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_FLOAT); 
     memcpy(PyArray_DATA(obj),&((*(arg8))[0]),sizeof(float)*length);	 
     delete arg8; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -3315,21 +3317,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg6))[0]),sizeof(int)*length);	 
     delete arg6; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg7)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg7))[0]),sizeof(int)*length);	 
     delete arg7; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg8)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_DOUBLE); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_DOUBLE); 
     memcpy(PyArray_DATA(obj),&((*(arg8))[0]),sizeof(double)*length);	 
     delete arg8; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -3580,21 +3582,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(int)*length);	 
     delete arg11; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg12)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg12))[0]),sizeof(int)*length);	 
     delete arg12; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg13)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_FLOAT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_FLOAT); 
     memcpy(PyArray_DATA(obj),&((*(arg13))[0]),sizeof(float)*length);	 
     delete arg13; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -3802,21 +3804,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(int)*length);	 
     delete arg11; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg12)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg12))[0]),sizeof(int)*length);	 
     delete arg12; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg13)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_DOUBLE); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_DOUBLE); 
     memcpy(PyArray_DATA(obj),&((*(arg13))[0]),sizeof(double)*length);	 
     delete arg13; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -4088,21 +4090,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg6))[0]),sizeof(int)*length);	 
     delete arg6; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg7)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg7))[0]),sizeof(int)*length);	 
     delete arg7; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg8)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_FLOAT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_FLOAT); 
     memcpy(PyArray_DATA(obj),&((*(arg8))[0]),sizeof(float)*length);	 
     delete arg8; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -4211,21 +4213,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg6)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg6))[0]),sizeof(int)*length);	 
     delete arg6; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg7)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg7))[0]),sizeof(int)*length);	 
     delete arg7; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg8)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_DOUBLE); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_DOUBLE); 
     memcpy(PyArray_DATA(obj),&((*(arg8))[0]),sizeof(double)*length);	 
     delete arg8; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -4449,21 +4451,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg9)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg9))[0]),sizeof(int)*length);	 
     delete arg9; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg10)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg10))[0]),sizeof(int)*length);	 
     delete arg10; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_FLOAT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_FLOAT); 
     memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(float)*length);	 
     delete arg11; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
@@ -4626,21 +4628,21 @@
   resultobj = SWIG_Py_Void();
   {
     int length = (arg9)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg9))[0]),sizeof(int)*length);	 
     delete arg9; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg10)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_INT); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_INT); 
     memcpy(PyArray_DATA(obj),&((*(arg10))[0]),sizeof(int)*length);	 
     delete arg10; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 
   }
   {
     int length = (arg11)->size(); 
-    PyObject *obj = PyArray_FromDims(1, &length, PyArray_DOUBLE); 
+    PyObject *obj = PyArray_FromDims(1, &length,PyArray_DOUBLE); 
     memcpy(PyArray_DATA(obj),&((*(arg11))[0]),sizeof(double)*length);	 
     delete arg11; 
     resultobj = helper_appendToTuple( resultobj, (PyObject *)obj ); 

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-08-21 22:12:51 UTC (rev 3252)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-08-22 15:55:13 UTC (rev 3253)
@@ -1,11 +1,12 @@
-from scipy import ones,zeros,rand,array,array_split,hstack,transpose,sum,ones_like,sqrt
-from scipy.sparse import spidentity
+from scipy import ones,zeros,rand,array,array_split,hstack,transpose,sum,ones_like,sqrt,concatenate
+from scipy.sparse import spidentity,csr_matrix,coo_matrix
 from numpy.linalg import norm
+from numpy import zeros_like,arange,inner,diff,ravel
 
 import pydec
 from pydec import diag_sparse,inf_norm, mls_polynomial_coeffs,polynomial_smoother
      
-from multigrid import sa_interpolation,rs_interpolation
+from multigrid import sa_interpolation,rs_interpolation,sa_constant_interpolation,sa_no_threshold
 import multigrid
 import multigridtools
 from relaxation import gauss_seidel,jacobi 
@@ -312,7 +313,6 @@
         
         self[-1].coarse_solver = coarse_grid_solver(A,self.opts.sub_options('coarse:'))
         self[-1].A = A
-
         if self.opts['smoother: type'] == 'jacobi':
             omegas = []
             for lvl in self:
@@ -324,4 +324,400 @@
             self.opts['smoother: omega'] = omegas
         
 
+class multilevel_solver2:
+    def __init__(self,As,Ps,options=None):
+        self.As = As
+        self.Ps = Ps
+        self.ops = options
+            
+    def solve(self,b, x0=None, tol=1e-5, maxiter=100, callback=None, return_residuals=False):
 
+        if x0 is None:
+            x = zeros(b.shape,max(self.A.dtype,b.dtype))
+        else:
+            x = array(x0)
+
+        self.__solve(0,x,b)
+
+        return x
+        
+    def __solve(self,lvl,x,b):
+
+        A = self.As[lvl]
+        
+        if len(self.As) == 1:
+            x[:] = scipy.linalg.solve(A.todense(),b)
+            return x
+
+
+        self.__smooth(lvl,x,b,which='pre')
+            
+        residual = b - A*x
+
+        coarse_x = zeros((self.As[lvl+1].shape[0]))
+        coarse_b = self.Ps[lvl].T * residual
+        
+        if lvl == len(self.As) - 2:
+            pass
+            coarse_x[:] = scipy.linalg.solve(self.As[-1].todense(),coarse_b)
+            #coarse_x[:] = self[-1].coarse_solver.solve(coarse_b)   #next level is coarsest 
+        else:   
+            self.__solve(lvl+1,coarse_x,coarse_b)
+                
+        x += self.Ps[lvl] * coarse_x
+
+        self.__smooth(lvl,x,b,which='post')
+
+
+    def __smooth(self,lvl,x,b,which):
+        A = self.As[lvl]
+        if which == 'pre':
+            gauss_seidel(A,x,b,iterations=1,sweep="forward")
+        else:
+            gauss_seidel(A,x,b,iterations=1,sweep="backward")
+
+
+def inf_norm(A):
+    return abs(A).sum(axis=1).max()  #max abs row sum
+
+def fit_candidate(I,x):
+    """
+    For each aggregate in I (i.e. each column of I) compute vector R and 
+    sparse matrix Q (having the sparsity of I) such that the following holds:
+
+        Q*R = x     and   Q^T*Q = I
+
+    In otherwords, find a prolongator Q with orthonormal columns so that
+    x is represented exactly on the coarser level by R.
+    """
+    Q = csr_matrix((x.copy(),I.indices,I.indptr),dims=I.shape,check=False)
+    R = sqrt(numpy.ravel(csr_matrix((x*x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0)))  #column 2-norms  
+    Q.data *= (1.0/R)[Q.indices]
+    print "norm(Q*R - x)",linalg.norm(Q*R - x)
+    return Q,R
+
+
+def scaled_columns_csr(A,scales):
+    scales = numpy.ravel(scales)
+    A = A.copy()
+    A.data *= scales[A.indices]
+    return A
+
+def orthonormalize_candidate(I,x,basis):
+    Px = csr_matrix((x,I.indices,I.indptr),dims=I.shape,check=False) 
+    Rs = []
+    #othogonalize columns of Px against other candidates 
+    for b in basis:
+        Pb = csr_matrix((b,I.indices,I.indptr),dims=I.shape,check=False)
+        R  = ravel(csr_matrix((Pb.data*Px.data,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0)) # columnwise projection of Px on Pb
+        Px.data -= R[I.indices] * Pb.data   #subtract component in b direction
+        Rs.append(R)
+
+    #filter columns here, set unused cols to 0, add to mask
+    
+    #normalize columns of Px
+    R = ravel(csr_matrix((x**x,I.indices,I.indptr),dims=I.shape,check=False).sum(axis=0))
+    Px.data *= (1.0/R)[I.indices]
+    Rs.append(R.reshape(-1,1))
+    return Rs
+
+def hstack_csr(A,B):
+    #OPTIMIZE THIS
+    assert(A.shape[0] == B.shape[0])
+    A = A.tocoo()
+    B = B.tocoo()
+    I = concatenate((A.row,B.row))
+    J = concatenate((A.col,B.col+A.shape[1]))
+    V = concatenate((A.data,B.data))
+    return coo_matrix((V,(I,J)),dims=(A.shape[0],A.shape[1]+B.shape[1])).tocsr()
+
+
+def vstack_csr(A,B):
+    #OPTIMIZE THIS
+    assert(A.shape[1] == B.shape[1])
+    A = A.tocoo()
+    B = B.tocoo()
+    I = concatenate((A.row,B.row+A.shape[0]))
+    J = concatenate((A.col,B.col))
+    V = concatenate((A.data,B.data))
+    return coo_matrix((V,(I,J)),dims=(A.shape[0]+B.shape[0],A.shape[1])).tocsr()
+    
+
+
+def orthonormalize_prolongator(P_l,x_l,W_l,W_m):
+    """
+     
+    """
+    X = csr_matrix((x_l,W_l.indices,W_l.indptr),dims=W_l.shape,check=False)  #candidate prolongator (assumes every value from x is used)
+    
+    R = (P_l.T.tocsr() * X)  # R has at most 1 nz per row
+    X = X - P_l*R            # othogonalize X against P_l
+    
+    #DROP REDUNDANT COLUMNS FROM P (AND R?) HERE (NULL OUT R ACCORDINGLY?)    
+    #REMOVE CORRESPONDING COLUMNS FROM W_l AND ROWS FROM A_m ALSO
+    W_l_new = W_l
+    W_m_new = W_m
+
+    #normalize surviving columns of X
+    col_norms = ravel(sqrt(csr_matrix((X.data*X.data,X.indices,X.indptr),dims=X.shape,check=False).sum(axis=0)))
+    print "zero cols",sum(col_norms == 0)
+    print "small cols",sum(col_norms < 1e-8)
+    Xcopy = X.copy()
+    X.data *= (1.0/col_norms)[X.indices]
+
+    P_l_new = hstack_csr(P_l,X)
+
+
+    #check orthonormality
+    print "norm(P.T*P - I) ",scipy.linalg.norm((P_l_new.T * P_l_new - scipy.sparse.spidentity(P_l_new.shape[1])).data)
+    #assert(scipy.linalg.norm((P_l_new.T * P_l_new - scipy.sparse.spidentity(P_l_new.shape[1])).data)<1e-8)
+    
+    x_m = zeros(P_l_new.shape[1],dtype=x_l.dtype)
+    x_m[:P_l.shape[1]][diff(R.indptr).astype('bool')] = R.data
+    x_m[P_l.shape[1]:] = col_norms
+
+    print "||x_l - P_l*x_m||",scipy.linalg.norm(P_l_new* x_m - x_l) #see if x_l is represented exactly
+
+    return P_l_new,x_m,W_l,W_m
+
+
+
+def prolongation_smoother(A):
+    omega = (4.0/3.0)/inf_norm(A)
+    S = (spidentity(A.shape[0]).T  - omega*A)
+    return S
+
+
+def smoothed_prolongator(P,A):
+    #just use Richardson for now
+    omega = 4.0/(3.0*inf_norm(A))
+    return P - omega*(A*P)
+ 
+
+
+
+def sa_hierarchy(A,Ws,x):
+    """
+    Construct multilevel hierarchy using Smoothed Aggregation
+        Inputs:
+          A  - matrix
+          Is - list of constant prolongators
+          x  - "candidate" basis function to be approximated
+        Ouputs:
+          (As,Is,Ps) - tuple of lists
+                  - As - [A, Ps[0].T*A*Ps[0], Ps[1].T*A*Ps[1], ... ]
+                  - Is - smoothed prolongators 
+                  - Ps - tentative prolongators
+    """
+    Ps = []
+    Is = []
+    As = [A]
+
+    for W in Ws:
+        P,x = fit_candidate(W,x)
+        I   = smoothed_prolongator(P,A)  
+        A   = I.T.tocsr() * A * I
+        As.append(A)
+        Ps.append(P)
+        Is.append(I)
+    return As,Is,Ps
+         
+def make_bridge(I,N):
+    tail = I.indptr[-1].repeat(N - I.shape[0])
+    ptr = concatenate((I.indptr,tail))
+    return csr_matrix((I.data,I.indices,ptr),dims=(N,I.shape[1]),check=False)
+
+class adaptive_sa_solver:
+    def __init__(self,A,options=None):
+        self.A = A
+
+        self.Rs = [] 
+        self.__construct_hierarchy(A)
+
+    def __construct_hierarchy(self,A):
+        #if self.A.shape[0] <= self.opts['coarse: max size']:
+        #    raise ValueError,'small matrices not handled yet'
+        
+        x,AggOps = self.__initialization_stage(A) #first candidate
+        Ws = AggOps
+        
+        #x[:] = 1 #TEST 
+
+        self.candidates = [x]
+        
+        #create SA using x here
+        As,Is,Ps = sa_hierarchy(A,Ws,x)
+
+        for i in range(0):
+            x           = self.__develop_candidate(A,As,Is,Ps,Ws,AggOps)
+            #x[:]  = arange(x.shape[0])
+            #x[x.shape[0]/2:] = 2*x[x.shape[0]/2] - x[x.shape[0]/2:]
+            As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)
+             
+            self.candidates.append(x)
+        
+        #As,Is,Ps = sa_hierarchy(A,AggOps,x)  #TESTING
+        self.Ps = Ps 
+        self.solver = multilevel_solver2(As,Is)
+                
+
+
+    def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps):
+        x = rand(A.shape[0])
+        b = zeros_like(x)
+
+        #x[:] = 1 #TEST
+
+        mu = 5
+ 
+        solver = multilevel_solver2(As,Is)
+
+        for n in range(mu):
+            x = solver.solve(b, x0=x, tol=1e-8, maxiter=1)
+        #TEST FOR CONVERGENCE HERE
+ 
+        A_l,P_l,W_l,x_l = As[0],Ps[0],Ws[0],x
+        
+        temp_Is = []
+        for i in range(len(As) - 2):
+            P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])    
+ 
+            I_l_new = smoothed_prolongator(P_l_new,A_l)
+            A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
+            bridge = make_bridge(Is[i+1],A_m_new.shape[0]) 
+ 
+            temp_solver = multilevel_solver2( [A_m_new] + As[i+2:], [bridge] + Is[i+2:] )
+ 
+            for n in range(mu):
+                x_m = temp_solver.solve(zeros_like(x_m), x0=x_m, tol=1e-8, maxiter=1)
+ 
+            temp_Is.append(I_l_new)
+ 
+            W_l = vstack_csr(Ws[i+1],W_m_new)  #prepare for next iteration
+            A_l = A_m_new
+            x_l = x_m
+            P_l = make_bridge(Ps[i+1],A_m_new.shape[0])
+ 
+        x = x_l
+        for I in reversed(temp_Is):
+            x = I*x
+        
+        return x
+           
+
+    def __augment_cycle(self,A,As,Ps,Ws,AggOps,x):
+        #As,Is,Ps,Ws = self.__augment_cycle(A,Ps,Ws,AggOps,x)
+
+        #make a new cycle using the new candidate
+        A_l,P_l,W_l,x_l = As[0],Ps[0],AggOps[0],x            
+        
+        new_As,new_Is,new_Ps,new_Ws = [A],[],[],[AggOps[0]]
+
+        for i in range(len(As) - 2):
+            P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, AggOps[i+1])
+
+            I_l_new = smoothed_prolongator(P_l_new,A_l)
+            A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
+            W_m_new = vstack_csr(Ws[i+1],W_m_new)
+
+            new_As.append(A_m_new)
+            new_Ws.append(W_m_new)
+            new_Is.append(I_l_new)
+            new_Ps.append(P_l_new)
+
+            #prepare for next iteration
+            W_l = W_m_new 
+            A_l = A_m_new
+            x_l = x_m
+            P_l = make_bridge(Ps[i+1],A_m_new.shape[0]) 
+        
+        P_l_new, x_m, W_l_new, W_m_new = orthonormalize_prolongator(P_l, x_l, W_l, csr_matrix((P_l.shape[1],1)))
+        I_l_new = smoothed_prolongator(P_l_new,A_l)
+        A_m_new = I_l_new.T.tocsr() * A_l * I_l_new
+
+        new_As.append(A_m_new)
+        new_Is.append(I_l_new)
+        new_Ps.append(P_l_new)
+
+        return new_As,new_Is,new_Ps,new_Ws
+
+
+    def __initialization_stage(self,A):
+        max_levels = 10
+        max_coarse = 50
+
+        AggOps = []
+        Is     = []
+
+        # aSA parameters 
+        mu      = 5    # number of test relaxation iterations
+        epsilon = 0.1  # minimum acceptable relaxation convergence factor 
+
+        #step 1
+        A_l = A
+        x   = scipy.rand(A_l.shape[0])
+        skip_f_to_i = False
+        
+        #step 2
+        b = zeros_like(x)
+        gauss_seidel(A_l,x,b,iterations=mu)
+        #step 3
+        #test convergence rate here 
+       
+        while len(AggOps) < max_levels  and A_l.shape[0] > max_coarse:
+            W_l   = sa_constant_interpolation(A_l)                        #step 4b 
+            #W_l   = sa_no_threshold(A_l)                        #step 4b  TEST
+            P_l,x = fit_candidate(W_l,x)                                  #step 4c  
+            I_l   = smoothed_prolongator(P_l,A_l)                         #step 4d
+            A_l   = I_l.T.tocsr() * A_l * I_l                             #step 4e
+            
+            AggOps.append(W_l)
+            Is.append(I_l)
+
+            if A_l.shape <= max_coarse:  break
+
+            if not skip_f_to_i:
+                print "."
+                x_hat = x.copy()                                        #step 4g
+                gauss_seidel(A_l,x,zeros_like(x),iterations=mu)         #step 4h
+                x_A_x = inner(x,A_l*x) 
+                if (x_A_x/inner(x_hat,A_l*x_hat))**(1.0/mu) < epsilon:  #step 4i
+                    print "sufficient convergence, skipping"
+                    skip_f_to_i = True
+                    if x_A_x == 0:       
+                        x = x_hat  #need to restore x
+        
+        #update fine-level candidate
+        for I in reversed(Is):
+            x = I * x
+                
+        #gauss_seidel(A,x,zeros_like(x),iterations=mu)    #TEST
+    
+        #x[:] = 1 #TEST
+    
+        return x,AggOps  #first candidate,aggregation
+
+
+
+from scipy import *
+from pydec import diag_sparse
+from multigrid import poisson_problem,poisson_problem1D
+#A = poisson_problem(100).T
+A = poisson_problem1D(100).T
+D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+A = D * A * D
+#A = A*A
+#A = io.mmread("nos2.mtx").tocsr()
+asa = adaptive_sa_solver(A)
+x = rand(A.shape[0])
+b = zeros_like(x)
+
+resid = []
+
+for n in range(50):
+    x = asa.solver.solve(b,x)
+    resid.append(linalg.norm(A*x))
+
+
+
+

Modified: trunk/scipy/sandbox/multigrid/simple_test.py
===================================================================
--- trunk/scipy/sandbox/multigrid/simple_test.py	2007-08-21 22:12:51 UTC (rev 3252)
+++ trunk/scipy/sandbox/multigrid/simple_test.py	2007-08-22 15:55:13 UTC (rev 3253)
@@ -2,7 +2,7 @@
 from multigrid import *
 from scipy import *
 
-A = poisson_problem(300).T
+A = poisson_problem(100).T
 s = scalar_solver(A)
 b = rand(A.shape[0])
 x,res = s.solve(b,return_residuals=True)


Property changes on: trunk/scipy/sandbox/multigrid/tests
___________________________________________________________________
Name: svn:ignore
   + *.so
*.bak
*.pyc
ignore.txt





More information about the Scipy-svn mailing list