# [Scipy-svn] r3432 - trunk/scipy/sandbox/multigrid

scipy-svn@scip... scipy-svn@scip...
Wed Oct 10 13:06:02 CDT 2007

```Author: wnbell
Date: 2007-10-10 13:06:00 -0500 (Wed, 10 Oct 2007)
New Revision: 3432

Modified:
trunk/scipy/sandbox/multigrid/utils.py
Log:

===================================================================
--- trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-10 17:12:37 UTC (rev 3431)
+++ trunk/scipy/sandbox/multigrid/adaptive.py	2007-10-10 18:06:00 UTC (rev 3432)
@@ -5,53 +5,9 @@
from relaxation import gauss_seidel
from multilevel import multilevel_solver
from sa import sa_constant_interpolation,sa_fit_candidates
-#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.
-##    """
-##    x = asarray(x)
-##    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]  #normalize columns of Q
-##
-##    #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])
@@ -327,7 +283,6 @@
for A_l,I in reversed(zip(As[1:],Is)):
gauss_seidel(A_l,x,zeros_like(x),iterations=mu,sweep='symmetric')         #TEST
x = I * x
-
gauss_seidel(A,x,b,iterations=mu)         #TEST

return x,AggOps  #first candidate,aggregation
@@ -344,28 +299,30 @@

-#A = A*A
+
#D = diag_sparse(1.0/sqrt(10**(12*rand(A.shape[0])-6))).tocsr()
#A = D * A * D
-#x = arange(A.shape[0]).astype('d') + 1
-scipy.random.seed(0)  #TEST
+
+
+scipy.random.seed(0)  #make tests repeatable
x = rand(A.shape[0])
-b = zeros_like(x)
+b = A*rand(A.shape[0])

print "solving"
-#x_sol,residuals = asa.solver.solve(b,x,tol=1e-8,maxiter=30,return_residuals=True)
-if True:
+if False:
x_sol,residuals = asa.solver.solve(b,x0=x,maxiter=10,tol=1e-12,return_residuals=True)
else:
residuals = []
residuals.append(linalg.norm(b - A*x))
A.psolve = asa.solver.psolve
+
residuals = array(residuals)/residuals[0]
+
print "residuals ",residuals
print "mean convergence factor",(residuals[-1]/residuals[0])**(1.0/len(residuals))
print "last convergence factor",residuals[-1]/residuals[-2]

Modified: trunk/scipy/sandbox/multigrid/utils.py
===================================================================
--- trunk/scipy/sandbox/multigrid/utils.py	2007-10-10 17:12:37 UTC (rev 3431)
+++ trunk/scipy/sandbox/multigrid/utils.py	2007-10-10 18:06:00 UTC (rev 3432)
@@ -1,7 +1,9 @@

-import numpy,scipy,scipy.sparse
+import numpy
+import scipy
from numpy import ravel,arange
+from scipy.linalg import norm
from scipy.sparse import isspmatrix,isspmatrix_csr,isspmatrix_csc, \
csr_matrix,csc_matrix,extract_diagonal

@@ -11,7 +13,7 @@
Approximate the spectral radius of a symmetric matrix using ARPACK
"""
from scipy.sandbox.arpack import eigen
-    return eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0]
+    return norm(eigen(A, k=1, ncv=10, which='LM', maxiter=maxiter, tol=tol, return_eigenvectors=False)[0])

```