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

scipy-svn@scip... scipy-svn@scip...
Fri Sep 7 14:56:01 CDT 2007


Author: wnbell
Date: 2007-09-07 14:55:57 -0500 (Fri, 07 Sep 2007)
New Revision: 3307

Added:
   trunk/scipy/sandbox/multigrid/adaptive.py
   trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
   trunk/scipy/sandbox/multigrid/tests/test_utils.py
Modified:
   trunk/scipy/sandbox/multigrid/coarsen.py
   trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
   trunk/scipy/sandbox/multigrid/multilevel.py
   trunk/scipy/sandbox/multigrid/relaxation.py
   trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
   trunk/scipy/sandbox/multigrid/utils.py
Log:
added test cases for utils and coarsening
added preliminary adaptive SA code



Added: trunk/scipy/sandbox/multigrid/adaptive.py
===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,371 @@
+import numpy,scipy,scipy.sparse
+from numpy import sqrt,ravel,diff,zeros,zeros_like,inner,concatenate
+from scipy.sparse import csr_matrix,coo_matrix
+
+from relaxation import gauss_seidel
+from multilevel import multilevel_solver
+from coarsen import sa_constant_interpolation
+from utils import infinity_norm
+
+
+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(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(R)",scipy.linalg.norm(R)
+    #print "min(R),max(R)",min(R),max(R)
+    #print "infinity_norm(Q.T*Q - I) ",infinity_norm((Q.T.tocsr() * Q - scipy.sparse.spidentity(Q.shape[1])))
+    #print "norm(Q*R - x)",scipy.linalg.norm(Q*R - x)
+    #print "norm(x - Q*Q.Tx)",scipy.linalg.norm(x - Q*(Q.T*x))
+    return Q,R
+
+
+##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 smoothed_prolongator(P,A):
+    #just use Richardson for now
+    #omega = 4.0/(3.0*infinity_norm(A))
+    #return P - omega*(A*P)
+    #return P
+    D = diag_sparse(A)
+    D_inv_A = diag_sparse(1.0/D)*A
+    omega = 4.0/(3.0*infinity_norm(D_inv_A))
+    D_inv_A *= omega
+    return P - D_inv_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
+    """
+    As = [A]
+    Is = []
+    Ps = []
+
+    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]
+        #self.candidates = [1.0/D.data]
+
+        #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)
+            #if i == 0:
+            #    x = arange(20).repeat(20).astype(float)
+            #elif i == 1:
+            #    x = arange(20).repeat(20).astype(float)
+            #    x = numpy.ravel(transpose(x.reshape((20,20))))
+
+            As,Is,Ps,Ws = self.__augment_cycle(A,As,Ps,Ws,AggOps,x)
+             
+            self.candidates.append(x)
+        
+        self.Ps = Ps 
+        self.solver = multilevel_solver(As,Is)
+        self.AggOps = AggOps
+                
+
+
+    def __develop_candidate(self,A,As,Is,Ps,Ws,AggOps):
+        x = scipy.rand(A.shape[0])
+        b = zeros_like(x)
+
+        
+        #x = arange(200).repeat(200).astype(float)
+        #x[:] = 1 #TEST
+
+        mu = 5
+ 
+        solver = multilevel_solver(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_solver( [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):
+        #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 
+        
+        scipy.random.seed(0) 
+
+        #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 
+      
+        As = [A]
+
+        while len(AggOps) + 1 < max_levels  and A_l.shape[0] > max_coarse:
+            W_l   = sa_constant_interpolation(A_l,epsilon=0.08*0.5**(len(AggOps)-1))           #step 4b  #TEST
+            #W_l   = sa_constant_interpolation(A_l,epsilon=0)           #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)
+            As.append(A_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 A_l,I in reversed(zip(As[1:],Is)):
+            gauss_seidel(A_l,x,zeros_like(x),iterations=mu)         #TEST
+            x = I * x
+        
+        gauss_seidel(A,x,b,iterations=mu)         #TEST
+    
+        return x,AggOps  #first candidate,aggregation
+
+
+
+
+from scipy import *
+from utils import diag_sparse
+from multilevel import poisson_problem1D,poisson_problem2D
+#A = poisson_problem2D(100)
+A = io.mmread("tests/sample_data/laplacian_40_3dcube.mtx").tocsr()
+
+#A = A*A
+#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
+#A = D * A * D
+#A = io.mmread("nos2.mtx").tocsr()
+asa = adaptive_sa_solver(A)
+x = rand(A.shape[0])
+b = zeros_like(x)
+
+
+print "solving"
+x_sol,residuals = asa.solver.solve(b,x,tol=1e-12,maxiter=30,return_residuals=True)
+residuals = array(residuals)/residuals[0]
+print "residuals ",residuals
+
+print asa.solver
+
+print "constant Rayleigh quotient",dot(ones(A.shape[0]),A*ones(A.shape[0]))/float(A.shape[0])
+
+for c in asa.candidates:
+    print "candidate Rayleigh quotient",dot(c,A*c)/dot(c,c)
+
+
+
+##W = asa.AggOps[0]*asa.AggOps[1]
+##pcolor((W * rand(W.shape[1])).reshape((200,200)))
+
+def plot2d(x):
+    from pylab import pcolor
+    pcolor(x.reshape(sqrt(len(x)),sqrt(len(x))))
+    

Modified: trunk/scipy/sandbox/multigrid/coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/coarsen.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/coarsen.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -1,10 +1,9 @@
-from scipy import *
 
 import multigridtools
 import scipy
 import numpy
     
-from utils import diag_sparse,inf_norm
+from utils import diag_sparse,infinity_norm
 
 
 def rs_strong_connections(A,theta):
@@ -33,37 +32,45 @@
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
 
     Sp,Sj,Sx = multigridtools.sa_strong_connections(A.shape[0],epsilon,A.indptr,A.indices,A.data)
+
     return scipy.sparse.csr_matrix((Sx,Sj,Sp),A.shape)
 
-def sa_constant_interpolation(A,epsilon=None):
+def sa_constant_interpolation(A,epsilon):
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
     
-    if epsilon is not None:
-        S = sa_strong_connections(A,epsilon)
-    else:
-        S = A
-    
+    S = sa_strong_connections(A,epsilon)
+
+    #S.ensure_sorted_indices()
+
     #tentative (non-smooth) interpolation operator I
-    Ij = multigridtools.sa_get_aggregates(A.shape[0],S.indptr,S.indices)
-    Ip = numpy.arange(len(Ij)+1)
-    Ix = numpy.ones(len(Ij))
+    Pj = multigridtools.sa_get_aggregates(S.shape[0],S.indptr,S.indices)
+    Pp = numpy.arange(len(Pj)+1)
+    Px = numpy.ones(len(Pj))
     
-    return scipy.sparse.csr_matrix((Ix,Ij,Ip))
+    return scipy.sparse.csr_matrix((Px,Pj,Pp))
 
-
+##def sa_smoother(A,S,omega):
+##    Bp,Bj,Bx = multigridtools.sa_smoother(A.shape[0],omega,A.indptr,A.indices,A.data,S.indptr,S.indices,S.data)
+##
+##    return csr_matrix((Bx,Bj,Bp),dims=A.shape)
+    
 def sa_interpolation(A,epsilon,omega=4.0/3.0):
     if not scipy.sparse.isspmatrix_csr(A): raise TypeError('expected sparse.csr_matrix')
+   
+    P  = sa_constant_interpolation(A,epsilon)
+
+##    As = sa_strong_connections(A,epsilon)
+##    S  = sa_smoother(A,S,omega)
     
-    I = sa_constant_interpolation(A,epsilon)
 
     D_inv = diag_sparse(1.0/diag_sparse(A))       
     
     D_inv_A  = D_inv * A
-    D_inv_A *= -omega/inf_norm(D_inv_A)
+    D_inv_A *= omega/infinity_norm(D_inv_A)
 
-    P = I + (D_inv_A*I)  #same as P=S*I, (faster?)
-        
-    return P
+    I = P - (D_inv_A*P)  #same as I=S*P, (faster?)
+           
+    return I
 
 
 

Modified: trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h
===================================================================
--- trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/multigridtools/smoothed_aggregation.h	2007-09-07 19:55:57 UTC (rev 3307)
@@ -5,8 +5,8 @@
 #include <vector>
 #include <iterator>
 #include <assert.h>
+#include <cmath>
 
-
 //#define DEBUG
 
 
@@ -20,24 +20,22 @@
   Sp->push_back(0);
 
   //compute diagonal values
-  std::vector<T> diags(n_row);
+  std::vector<T> diags(n_row,T(0));
   for(int i = 0; i < n_row; i++){
     int row_start = Ap[i];
     int row_end   = Ap[i+1];
     for(int jj = row_start; jj < row_end; jj++){
-      if(Aj[jj] == i){
-	diags[i] = Ax[jj];
-	break;
-      }
+        if(Aj[jj] == i){
+        	diags[i] = Ax[jj];
+        	break;
+        }
     }    
   }
 
 #ifdef DEBUG
   for(int i = 0; i < n_row; i++){ assert(diags[i] > 0); }
 #endif
-    
 
-
   for(int i = 0; i < n_row; i++){
     int row_start = Ap[i];
     int row_end   = Ap[i+1];
@@ -45,14 +43,15 @@
     T eps_Aii = epsilon*epsilon*diags[i];
 
     for(int jj = row_start; jj < row_end; jj++){
-      const int&   j = Aj[jj];
-      const T&   Aij = Ax[jj];
+      const int   j = Aj[jj];
+      const T   Aij = Ax[jj];
 
       if(i == j){continue;}
 
-      if(Aij*Aij >= eps_Aii * diags[j]){
-	Sj->push_back(j);
-	Sx->push_back(Aij);
+      //  |A(i,j)| < epsilon * sqrt(|A(i,i)|*|A(j,j)|) 
+      if(Aij*Aij >= std::abs(eps_Aii * diags[j])){    
+          Sj->push_back(j);
+          Sx->push_back(Aij);
       }
     }
     Sp->push_back(Sj->size());
@@ -61,9 +60,9 @@
 
 
 void sa_get_aggregates(const int n_row,
-		       const int Ap[], const int Aj[],
-		       std::vector<int> * Bj){
-
+        		       const int Ap[], const int Aj[],
+		               std::vector<int> * Bj)
+{
   std::vector<int> aggregates(n_row,-1);
 
   int num_aggregates = 0;
@@ -72,21 +71,19 @@
   for(int i = 0; i < n_row; i++){
     if(aggregates[i] >= 0){ continue; } //already marked
 
-    const int& row_start = Ap[i];
-    const int& row_end   = Ap[i+1];
+    const int row_start = Ap[i];
+    const int row_end   = Ap[i+1];
     
-
     //Determine whether all neighbors of this node are free (not already aggregates)
     bool free_neighborhood = true;
     for(int jj = row_start; jj < row_end; jj++){
       if(aggregates[Aj[jj]] >= 0){
-	free_neighborhood = false;
-	break;
+    	free_neighborhood = false;
+	    break;
       }
     }    
     if(!free_neighborhood){ continue; } //bail out
 
-
     //Make an aggregate out of this node and its strong neigbors
     aggregates[i] = num_aggregates;
     for(int jj = row_start; jj < row_end; jj++){
@@ -96,52 +93,49 @@
   }
 
 
-
   //Pass #2
   std::vector<int> aggregates_copy(aggregates);
   for(int i = 0; i < n_row; i++){
     if(aggregates[i] >= 0){ continue; } //already marked
 
-    const int& row_start = Ap[i];
-    const int& row_end   = Ap[i+1];
+    const int row_start = Ap[i];
+    const int row_end   = Ap[i+1];
     
     for(int jj = row_start; jj < row_end; jj++){
-      const int& j = Aj[jj];
+      const int j = Aj[jj];
 
       if(aggregates_copy[j] >= 0){
-	aggregates[i] = aggregates_copy[j];
-	break;
+    	aggregates[i] = aggregates_copy[j];
+	    break;
       }
     }    
   }
 
 
-
   //Pass #3
   for(int i = 0; i < n_row; i++){
     if(aggregates[i] >= 0){ continue; } //already marked
-
-    const int& row_start = Ap[i];
-    const int& row_end   = Ap[i+1];
+  
+    const int row_start = Ap[i];
+    const int row_end   = Ap[i+1];
     
     aggregates[i] = num_aggregates;
 
     for(int jj = row_start; jj < row_end; jj++){
-      const int& j = Aj[jj];
+      const int j = Aj[jj];
 
       if(aggregates[j] < 0){ //unmarked neighbors
-	aggregates[j] = num_aggregates;
+    	aggregates[j] = num_aggregates;
       }
     }  
     num_aggregates++;
   }
 
-
 #ifdef DEBUG
   for(int i = 0; i < n_row; i++){ assert(aggregates[i] >= 0 && aggregates[i] < num_aggregates); }
 #endif
   
-  *Bj = aggregates;  
+  aggregates.swap(*Bj);  
 }
 
 

Modified: trunk/scipy/sandbox/multigrid/multilevel.py
===================================================================
--- trunk/scipy/sandbox/multigrid/multilevel.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/multilevel.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -10,9 +10,9 @@
 
 from coarsen import sa_interpolation,rs_interpolation
 from relaxation import gauss_seidel,jacobi 
+from utils import infinity_norm
 
 
-
 def poisson_problem1D(N):
     """
     Return a sparse CSR matrix for the 1d poisson problem
@@ -21,7 +21,7 @@
     """
     D = 2*numpy.ones(N)
     O =  -numpy.ones(N)
-    return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocsr()
+    return scipy.sparse.spdiags([D,O,O],[0,-1,1],N,N).tocoo().tocsr() #eliminate zeros
 
 def poisson_problem2D(N):
     """
@@ -33,7 +33,8 @@
     T =  -numpy.ones(N*N)
     O =  -numpy.ones(N*N)
     T[N-1::N] = 0
-    return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocsr()
+    return scipy.sparse.spdiags([D,O,T,T,O],[0,-N,-1,1,N],N*N,N*N).tocoo().tocsr() #eliminate zeros
+    
 
 def ruge_stuben_solver(A,max_levels=10,max_coarse=500):
     """
@@ -41,8 +42,8 @@
     
         References:
             "Multigrid"
-            Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
-            See Appendix A
+                Trottenberg, U., C. W. Oosterlee, and Anton Schuller. San Diego: Academic Press, 2001.
+                See Appendix A
     
     """
     As = [A]
@@ -58,7 +59,7 @@
         
     return multilevel_solver(As,Ps)
 
-def smoothed_aggregation_solver(A,max_levels=10,max_coarse=500):
+def smoothed_aggregation_solver(A,max_levels=10,max_coarse=500,epsilon=0.08):
     """
     Create a multilevel solver using Smoothed Aggregation (SA)
 
@@ -72,8 +73,9 @@
     Ps = []
     
     while len(As) < max_levels  and A.shape[0] > max_coarse:
-        P = sa_interpolation(A,epsilon=0.08*0.5**(len(As)-1))
-        
+        P = sa_interpolation(A,epsilon=epsilon*0.5**(len(As)-1))
+        #P = sa_interpolation(A,epsilon=0.0)
+
         A = (P.T.tocsr() * A) * P     #galerkin operator
 
         As.append(A)
@@ -108,7 +110,10 @@
         """number of unknowns on all levels / number of unknowns on the finest level"""
         return sum([A.shape[0] for A in self.As])/float(self.As[0].shape[0])
     
-            
+
+    def psolve(self, b):
+        return self.solve(b,maxiter=1)
+
     def solve(self, b, x0=None, tol=1e-5, maxiter=100, callback=None, return_residuals=False):
         """
         TODO
@@ -122,12 +127,12 @@
 
         #TODO change use of tol (relative tolerance) to agree with other iterative solvers
         A = self.As[0]
-        residuals = [norm(b-A*x,2)]
+        residuals = [scipy.linalg.norm(b-A*x)]
 
         while len(residuals) <= maxiter and residuals[-1]/residuals[0] > tol:
             self.__solve(0,x,b)
 
-            residuals.append(scipy.linalg.norm(b-A*x,2))
+            residuals.append(scipy.linalg.norm(b-A*x))
 
             if callback is not None:
                 callback(x)
@@ -142,11 +147,11 @@
         A = self.As[lvl]
         
         if len(self.As) == 1:
-            x[:] = scipy.linalg.solve(A.todense(),b)
-            return x
+            x[:] = scipy.linsolve.spsolve(A,b)
+            return 
 
         self.presmoother(A,x,b)
-            
+
         residual = b - A*x                                
 
         coarse_x = zeros((self.As[lvl+1].shape[0]))
@@ -154,7 +159,9 @@
         
         if lvl == len(self.As) - 2:
             #direct solver on coarsest level
-            coarse_x[:] = scipy.linalg.solve(self.As[-1].todense(),coarse_b)
+            coarse_x[:] = scipy.linsolve.spsolve(self.As[-1],coarse_b)
+            #coarse_x[:] = scipy.linalg.cg(self.As[-1],coarse_b,tol=1e-12)[0]
+            #print "coarse residual norm",scipy.linalg.norm(coarse_b - self.As[-1]*coarse_x)
         else:   
             self.__solve(lvl+1,coarse_x,coarse_b)
                 
@@ -165,26 +172,32 @@
 
     def presmoother(self,A,x,b):
         gauss_seidel(A,x,b,iterations=1,sweep="forward")
+        #x += 4.0/(3.0*infinity_norm(A))*(b - A*x)
     
     def postsmoother(self,A,x,b):
-        gauss_seidel(A,x,b,iterations=1,sweep="backward")
+        gauss_seidel(A,x,b,iterations=1,sweep="forward")
+        #gauss_seidel(A,x,b,iterations=1,sweep="backward")
+        #x += 4.0/(3.0*infinity_norm(A))*(b - A*x)
 
 
 
 if __name__ == '__main__':
     from scipy import *
     A = poisson_problem2D(200)
+    #A = io.mmread("rocker_arm_surface.mtx").tocsr()
+
     ml = smoothed_aggregation_solver(A)
     #ml = ruge_stuben_solver(A)
+
     x = rand(A.shape[0])
     b = zeros_like(x)
+    #b = rand(A.shape[0])
     
-    resid = []
-    
-    for n in range(10):
-        x = ml.solve(b,x,maxiter=1)
-        resid.append(linalg.norm(A*x))
+    x_sol,residuals = ml.solve(b,x0=x,maxiter=40,tol=1e-10,return_residuals=True)
 
+    residuals = array(residuals)/residuals[0]
 
+    print residuals
 
 
+

Modified: trunk/scipy/sandbox/multigrid/relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/relaxation.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/relaxation.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -13,6 +13,12 @@
          iterations - number of iterations to perform (default: 1)
          sweep      - slice of unknowns to relax (default: all in forward direction)
     """ 
+    if A.shape[0] != A.shape[1]:
+        raise ValueError,'expected symmetric matrix'
+
+    if A.shape[1] != len(x) or len(x) != len(b):
+        raise ValueError,'unexpected number of unknowns'
+
     if sweep == 'forward':
         row_start,row_stop,row_step = 0,len(x),1
     elif sweep == 'backward':

Added: trunk/scipy/sandbox/multigrid/tests/test_coarsen.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_coarsen.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_coarsen.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,160 @@
+from numpy.testing import *
+
+from numpy import sqrt,empty,ones,arange,array_split
+from scipy import rand
+from scipy.sparse import spdiags,csr_matrix,lil_matrix
+import numpy
+
+set_package_path()
+import scipy.multigrid
+from scipy.multigrid.coarsen import sa_strong_connections,sa_constant_interpolation
+from scipy.multigrid.multilevel import poisson_problem1D,poisson_problem2D
+restore_path()
+
+
+def reference_sa_strong_connections(A,epsilon):
+    A_coo = A.tocoo()
+    S = lil_matrix(A.shape)
+
+    for (i,j,v) in zip(A_coo.row,A_coo.col,A_coo.data):
+        if i == j: continue #skip diagonal
+
+        if abs(A[i,j]) >= epsilon*sqrt(abs(A[i,i])*abs(A[j,j])):
+            S[i,j] = v
+
+    return S.tocsr()
+
+
+# note that this method only tests the current implementation, not
+# all possible implementations
+def reference_sa_constant_interpolation(A,epsilon):
+    S = sa_strong_connections(A,epsilon)
+    S = array_split(S.indices,S.indptr[1:-1])
+
+    n = A.shape[0]
+
+    R = set(range(n))
+    j = 0
+
+    aggregates = empty(n,dtype=A.indices.dtype)
+    aggregates[:] = -1
+
+
+    # Pass #1
+    for i,row in enumerate(S):
+        Ni = set(row) | set([i])
+
+        if Ni.issubset(R):
+            R -= Ni
+            for x in Ni:
+                aggregates[x] = j
+            j += 1
+
+    # Pass #2
+    Old_R = R.copy()
+    for i,row in enumerate(S):
+        if i not in R: continue
+        
+        for x in row:
+            if x not in Old_R:
+                aggregates[i] = aggregates[x]
+                R.remove(i)
+                break
+
+
+    # Pass #3
+    for i,row in enumerate(S):
+        if i not in R: continue
+        Ni = set(row) | set([i])
+
+        for x in Ni:
+            if x in R:
+                aggregates[x] = j
+            j += 1
+
+    assert(len(R) == 0)
+
+    Pj = aggregates
+    Pp = arange(n+1)
+    Px = ones(n)
+ 
+    return csr_matrix((Px,Pj,Pp))
+
+class test_sa_strong_connections(NumpyTestCase):
+    def check_simple(self):
+        N = 4
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
+        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+       
+        N = 100
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).tocsr()
+        S = spdiags([ -ones(N),-ones(N)],[-1,1],N,N).tocsr()
+        assert_array_equal(sa_strong_connections(A,0.50).todense(),S.todense())   #all connections are strong
+        assert_array_equal(sa_strong_connections(A,0.51).todense(),0*S.todense()) #no connections are strong
+
+    def check_random(self):
+        numpy.random.seed(0)
+
+        for N in [2,3,5,10]:
+            A = csr_matrix(rand(N,N))
+            for epsilon in [0.0,0.1,0.5,0.8,1.0,10.0]:
+                S_result = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson1D(self):
+        for N in [2,3,5,7,10,11,19]:
+            A = poisson_problem1D(N)
+            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+                S_result   = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson2D(self):
+        for N in [2,3,5,7,10,11,19]:
+            A = poisson_problem2D(N)
+            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+                S_result   = sa_strong_connections(A,epsilon)
+                S_expected = reference_sa_strong_connections(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+##    def check_sample_data(self):
+##        for filename in all_matrices:
+##            A = open_matrix(filename)
+
+
+S_result = None
+S_expected = None
+class test_sa_constant_interpolation(NumpyTestCase):
+    def check_random(self):
+        numpy.random.seed(0)
+
+        for N in [2,3,5,10]:
+            A = csr_matrix(rand(N,N))
+            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson1D(self):
+        for N in [2,3,5,7,10,11,20,21,29,30]:
+            A = poisson_problem1D(N)
+            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+    def check_poisson2D(self):
+        for N in [2,3,5,7,10,11,20,21,29,30]:
+            A = poisson_problem2D(N)
+            for epsilon in [0.0,0.1,0.5,0.8,1.0]:
+                S_result   = sa_constant_interpolation(A,epsilon)
+                S_expected = reference_sa_constant_interpolation(A,epsilon)
+                assert_array_equal(S_result.todense(),S_expected.todense())
+
+
+if __name__ == '__main__':
+    NumpyTest().run()
+

Modified: trunk/scipy/sandbox/multigrid/tests/test_relaxation.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_relaxation.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_relaxation.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -7,6 +7,7 @@
 
 
 set_package_path()
+import scipy.multigrid
 from scipy.multigrid.relaxation import polynomial_smoother,gauss_seidel,jacobi
 restore_path()
 
@@ -36,7 +37,51 @@
         polynomial_smoother(A,x,b,[-0.14285714,  1., -2.])
         assert_almost_equal(x,x0 - 0.14285714*A*A*r + A*r - 2*r)
 
+    def check_jacobi(self):
+        N = 1
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        jacobi(A,x,b)
+        assert_almost_equal(x,array([0]))
 
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = zeros(N)
+        b = arange(N).astype(numpy.float64)
+        jacobi(A,x,b)
+        assert_almost_equal(x,array([0.0,0.5,1.0]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = zeros(N)
+        jacobi(A,x,b)
+        assert_almost_equal(x,array([0.5,1.0,0.5]))
+
+        N = 1
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = array([10])
+        jacobi(A,x,b)
+        assert_almost_equal(x,array([5]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        b = array([10,20,30])
+        jacobi(A,x,b)
+        assert_almost_equal(x,array([5.5,11.0,15.5]))
+
+        N = 3
+        A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T
+        x = arange(N).astype(numpy.float64)
+        x_copy = x.copy()
+        b = array([10,20,30])
+        jacobi(A,x,b,omega=1.0/3.0)
+        assert_almost_equal(x,2.0/3.0*x_copy + 1.0/3.0*array([5.5,11.0,15.5]))
+
+
     def check_gauss_seidel(self):
         N = 1
         A = spdiags([2*ones(N),-ones(N),-ones(N)],[0,-1,1],N,N).T

Added: trunk/scipy/sandbox/multigrid/tests/test_utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/tests/test_utils.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -0,0 +1,60 @@
+from numpy.testing import *
+
+import numpy
+import scipy
+from scipy import matrix,array,diag
+from scipy.sparse import csr_matrix
+
+
+set_package_path()
+from scipy.multigrid.utils import infinity_norm,diag_sparse
+restore_path()
+
+
+class test_utils(NumpyTestCase):
+    def check_infinity_norm(self):
+        A = matrix([[-4]])
+        assert_equal(infinity_norm(csr_matrix(A)),4)
+        
+        A = matrix([[1,0,-5],[-2,5,0]])
+        assert_equal(infinity_norm(csr_matrix(A)),7)
+
+        A = matrix([[0,1],[0,-5]])
+        assert_equal(infinity_norm(csr_matrix(A)),5)
+
+        A = matrix([[1.3,-4.7,0],[-2.23,5.5,0],[9,0,-2]])
+        assert_equal(infinity_norm(csr_matrix(A)),11)
+
+    def check_diag_sparse(self):
+        #check sparse -> array
+        A = matrix([[-4]])
+        assert_equal(diag_sparse(csr_matrix(A)),[-4])
+        
+        A = matrix([[1,0,-5],[-2,5,0]])
+        assert_equal(diag_sparse(csr_matrix(A)),[1,5])
+
+        A = matrix([[0,1],[0,-5]])
+        assert_equal(diag_sparse(csr_matrix(A)),[0,-5])
+
+        A = matrix([[1.3,-4.7,0],[-2.23,5.5,0],[9,0,-2]])
+        assert_equal(diag_sparse(csr_matrix(A)),[1.3,5.5,-2])
+
+        #check array -> sparse
+        A = matrix([[-4]])
+        assert_equal(diag_sparse(array([-4])).todense(),csr_matrix(A).todense())
+
+        A = matrix([[1,0],[0,5]])
+        assert_equal(diag_sparse(array([1,5])).todense(),csr_matrix(A).todense())
+
+        A = matrix([[0,0],[0,-5]])
+        assert_equal(diag_sparse(array([0,-5])).todense(),csr_matrix(A).todense())
+
+        A = matrix([[1.3,0,0],[0,5.5,0],[0,0,-2]])
+        assert_equal(diag_sparse(array([1.3,5.5,-2])).todense(),csr_matrix(A).todense())
+
+
+
+
+if __name__ == '__main__':
+    NumpyTest().run()
+

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-09-06 08:26:39 UTC (rev 3306)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-09-07 19:55:57 UTC (rev 3307)
@@ -6,17 +6,18 @@
                         csr_matrix,csc_matrix,extract_diagonal
 
 
-def inf_norm(A):
+def infinity_norm(A):
     """
     Infinity norm of a sparse matrix (maximum absolute row sum).  This serves 
     as an upper bound on spectral radius.
     """
     
-    if not isspmatrix_csr(A):
-        return ValueError,'expected csr_matrix'
-    
-    abs_A = csr_matrix((abs(A.data),A.indices,A.indptr),dims=A.shape,check=False)
-    return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+    if isspmatrix_csr(A) or isspmatrix_csc(A):
+        #avoid copying index and ptr arrays
+        abs_A = A.__class__((abs(A.data),A.indices,A.indptr),dims=A.shape,check=False)
+        return (abs_A * numpy.ones(A.shape[1],dtype=A.dtype)).max()
+    else:
+        return (abs(A) * numpy.ones(A.shape[1],dtype=A.dtype)).max()
 
 def diag_sparse(A):
     """



More information about the Scipy-svn mailing list