← Back to branch summary

~ubuntu-branches/ubuntu/precise/deal.ii/precise

« back to all changes in this revision

Viewing changes to contrib/umfpack/UMFPACK/MATLAB/lu_normest.m

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-05-08 23:13:50 UTC
  • Revision ID: james.westby@ubuntu.com-20090508231350-rrh1ltgi0tifabwc
Tags: upstream-6.2.0
Import upstream version 6.2.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/umfpack/UMFPACK/MATLAB/lu_normest.m
 
1
function rho = lu_normest (A, L, U)
 
2
%LU_NORMEST estimates norm (L*U-A, 1) without forming L*U-A
 
3
%
 
4
% Example:
 
5
%
 
6
%       rho = lu_normest (A, L, U)
 
7
%
 
8
% which estimates the computation of the 1-norm:
 
9
%
 
10
%       rho = norm (A-L*U, 1)
 
11
%
 
12
% Authors:  William W. Hager, Math Dept., Univ. of Florida
 
13
%       Timothy A. Davis, CISE Dept., Univ. of Florida
 
14
%       Gainesville, FL, 32611, USA.
 
15
%       based on normest1, contributed on November, 1997
 
16
%
 
17
% This code can be quite easily adapted to estimate the 1-norm of any
 
18
% matrix E, where E itself is dense or not explicitly represented, but the
 
19
% computation of E (and E') times a vector is easy.  In this case, our matrix
 
20
% of interest is:
 
21
%
 
22
%       E = A-L*U
 
23
%
 
24
% That is, L*U is the LU factorization of A, where A, L and U
 
25
% are sparse.  This code works for dense matrices A and L too,
 
26
% but it would not be needed in that case, since E is easy to compute
 
27
% explicitly.  For sparse A, L, and U, computing E explicitly would be quite
 
28
% expensive, and thus normest (A-L*U) would be prohibitive.
 
29
%
 
30
% For a detailed description, see Davis, T. A. and Hager, W. W.,
 
31
% Modifying a sparse Cholesky factorization, SIAM J. Matrix Analysis and
 
32
% Applications, 1999, vol. 20, no. 3, 606-627.
 
33
%
 
34
% See also normest
 
35
 
 
36
% The three places that the matrix-vector multiply E*x is used are highlighted.
 
37
% Note that E is never formed explicity.
 
38
 
 
39
% Copyright 1995-2007 by William W. Hager and Timothy A. Davis
 
40
 
 
41
[m n] = size (A) ;
 
42
 
 
43
if (m ~= n)
 
44
    % pad A, L, and U with zeros so that they are all square
 
45
    if (m < n)
 
46
        U = [ U ; (sparse (n-m,n)) ] ;
 
47
        L = [ L , (sparse (m,n-m)) ; (sparse (n-m,n)) ] ;
 
48
        A = [ A ; (sparse (n-m,n)) ] ;
 
49
    else
 
50
        U = [ U , (sparse (n,m-n)) ; (sparse (m-n,m)) ] ;
 
51
        L = [ L , (sparse (m,m-n)) ] ;
 
52
        A = [ A , (sparse (m,m-n)) ] ;
 
53
    end
 
54
end
 
55
 
 
56
[m n] = size (A) ;          %#ok
 
57
 
 
58
notvisited = ones (m, 1) ;  % nonvisited(j) is zero if j is visited, 1 otherwise
 
59
rho = 0 ;    % the global rho
 
60
 
 
61
for trial = 1:3 % {
 
62
 
 
63
   x = notvisited ./ sum (notvisited) ;
 
64
   rho1 = 0 ;    % the current rho for this trial
 
65
 
 
66
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
67
   %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
68
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
69
   Ex1 = (A*x) - L*(U*x) ;
 
70
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
71
 
 
72
   rho2 = norm (Ex1, 1) ;
 
73
 
 
74
   while rho2 > rho1 % {
 
75
 
 
76
        rho1 = rho2 ;
 
77
        y = 2*(Ex1 >= 0) - 1 ;
 
78
 
 
79
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
80
        %%% COMPUTE z = E'*y EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
81
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
82
        z = (A'*y) - U'*(L'*y) ;
 
83
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
84
 
 
85
        [zj, j] = max (abs (z .* notvisited)) ;
 
86
        j = j (1) ;
 
87
        if (abs (z (j)) > z'*x) % {
 
88
            x = zeros (m, 1) ;
 
89
            x (j) = 1 ;
 
90
            notvisited (j) = 0 ;
 
91
        else % } {
 
92
            break ;
 
93
        end % }
 
94
 
 
95
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
96
        %%% COMPUTE Ex1 = E*x EFFICIENTLY: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
97
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
98
        Ex1 = (A*x) - L*(U*x) ;
 
99
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
100
 
 
101
        rho2 = norm (Ex1, 1) ;
 
102
 
 
103
    end % }
 
104
 
 
105
    rho = max (rho, rho1) ;
 
106
 
 
107
end % }