1function B = spqr_null_mult (N,A,method)
2%SPQR_NULL_MULT multiplies a matrix by numerical null space from spqr_rank methods
3%
4%  Multiplies the matrix A by N, the orthonormal basis for
5%  a numerical null space produced by spqr_basic, spqr_null, spqr_pinv,
6%  or spqr_cod.  N can be stored either implicitly (see the above
7%  routines for details) or as an explicit matrix.
8%
9% Example of use:  B = spqr_null_mult(N,A,method) ;
10%    method = 0: N'*A  (the default)
11%    method = 1: N*A
12%    method = 2: A*N'
13%    method = 3: A*N
14% Also if N is stored implicitly then to create an explicit representation
15% of the orthonormal null space basis
16%    Nexplicit = spqr_null_mult(N,eye(size(N.X,2)),1) ;
17% creates a full matrix and
18%    Nexplicit = spqr_null_mult(N,speye(size(N.X,2)),1) ;
19% creates a sparse matrix.
20%
21% If N is stored implicitly then N can have fields:
22%    N.P contains a permutation vector (N.P is not always present)
23%    N.Q contains sparse Householder transforms from spqr
24%    N.X is a sparse matrix with orthonormal columns
25%
26% N.kind determines how the implicit orthonormal basis is represented.
27%   N.kind = 'P*Q*X':  the basis is N.P*N.Q*N.X
28%   N.kind = 'Q*P*X':  the basis is N.Q*N.P*N.X
29%   N.kind = 'Q*X':    the basis is N.Q*N.X, and N.P does not appear in N
30%
31% Examples
32%
33%   N = spqr_null (A) ;              % find nullspace of A
34%   AN = spqr_null_mult (N,A,3) ;    % compute A*N, which will have a low norm.
35%
36% See also spqr, spqr_basic, spqr_cod, spqr_null, spqr_pinv, spqr_ssi, spqr_ssp.
37
38% Copyright 2012, Leslie Foster and Timothy A Davis.
39
40% If N is stored implicitly, to potentially improve efficiency, code is
41% selected based on the number of rows and columns in A and N.
42
43if nargin < 3
44    method = 0 ;
45end
46[m n] = size (A) ;
47
48if isstruct (N)
49
50    %---------------------------------------------------------------------------
51    % N held implicitly
52    %---------------------------------------------------------------------------
53
54    p = size (N.X, 2) ;
55
56    switch method
57
58        case 0
59
60            %-------------------------------------------------------------------
61            % B = N'*A
62            %-------------------------------------------------------------------
63
64            switch N.kind
65
66                case 'Q*X'
67
68                    if n <= p
69                        % B = X' * ( Q' * A )
70                        B = spqr_qmult (N.Q, A, 0) ;
71                        B = N.X'*B ;
72                    else
73                        % B = ( X' * Q' ) * A
74                        B = spqr_qmult (N.Q, (N.X)', 2) ;
75                        B = B*A ;
76                    end
77
78                case 'Q*P*X'
79
80                    if n <= p
81                        % B = X' * P' * ( Q' * A ) ;
82                        B = spqr_qmult (N.Q,A,0) ;
83                        p (N.P) = 1:length (N.P) ;
84                        B = B(p,:) ;
85                        B = N.X'*B ;
86                    else
87                        % B = ( ( P * X )' * Q' ) * A ;
88                        B = N.X(N.P,:)';
89                        B = spqr_qmult (N.Q,B,2) ;
90                        B = B * A ;
91                    end
92
93                case 'P*Q*X'
94
95                    if n <= p
96                        % B = X' * ( Q' * ( P' * A ) ) ;
97                        p( N.P) = 1:length(N.P) ;
98                        B = A(p,:) ;
99                        B = spqr_qmult (N.Q,B,0) ;
100                        B = N.X'*B ;
101                    else
102                        % B = ( ( X' * Q' ) * P' ) * A ;
103                        B = spqr_qmult (N.Q, (N.X)', 2) ;
104                        B = B(:,N.P) ;
105                        B = B * A ;
106                    end
107
108                otherwise
109                    error ('unrecognized N struct') ;
110            end
111
112        case 1
113
114            %-------------------------------------------------------------------
115            % B = N*A
116            %-------------------------------------------------------------------
117
118            switch N.kind
119
120                case 'Q*X'
121
122                    % B = Q * ( X * A ) ;
123                    B = N.X * A ;
124                    B = spqr_qmult (N.Q,B,1) ;
125
126                case 'Q*P*X'
127
128                    % B = Q * ( P * ( X * A ) ) ;
129                    B = N.X * A ;
130                    B = B(N.P,:) ;
131                    B = spqr_qmult (N.Q,B,1) ;
132
133                case 'P*Q*X'
134
135                    % B = P * ( Q * ( X * A ) ) ;
136                    B = N.X * A ;
137                    B = spqr_qmult (N.Q,B,1) ;
138                    B = B(N.P,:) ;
139
140                otherwise
141                    error ('unrecognized N struct') ;
142            end
143
144        case 2
145
146            %-------------------------------------------------------------------
147            % B = A*N'
148            %-------------------------------------------------------------------
149
150            switch N.kind
151
152                case 'Q*X'
153
154                    % B = (A * X') * Q'
155                    B = A * (N.X)' ;
156                    B = spqr_qmult (N.Q,B,2) ;
157
158                case 'Q*P*X'
159
160                    % B = ( ( A * X' ) * P' ) * Q'
161                    B = A * (N.X)' ;
162                    B = B(:,N.P) ;
163                    B = spqr_qmult (N.Q,B,2) ;
164
165                case 'P*Q*X'
166
167                    % B = ( ( A * X' ) * Q' ) * P'
168                    B = A * (N.X)' ;
169                    B = spqr_qmult (N.Q,B,2) ;
170                    B = B(:,N.P) ;
171
172                otherwise
173                    error ('unrecognized N struct') ;
174            end
175
176        case 3
177
178            %-------------------------------------------------------------------
179            % B = A*N
180            %-------------------------------------------------------------------
181
182            switch N.kind
183
184                case 'Q*X'
185
186                    if m <= p
187                        % B = ( A * Q ) * X
188                        B = spqr_qmult (N.Q,A,3) ;
189                        B = B*N.X ;
190                    else
191                        % B = A * ( Q * X )
192                        B = spqr_qmult (N.Q, N.X, 1 ) ;
193                        B = A * B ;
194                    end
195
196                case 'Q*P*X'
197
198                    if m <= p
199                        % B = ( ( A * Q ) * P ) * X
200                        B = spqr_qmult (N.Q,A,3) ;
201                        p(N.P) = 1:length(N.P) ;
202                        B = B(:,p) ;
203                        B = B*N.X ;
204                    else
205                        % B = A * ( Q * ( P * X ) )
206                        B = N.X(N.P,:) ;
207                        B = spqr_qmult (N.Q, B, 1 ) ;
208                        B = A * B ;
209                    end
210
211                case 'P*Q*X'
212
213                    if m <= p
214                        % B = ( ( A * P ) * Q ) * X
215                        p(N.P) = 1:length(N.P) ;
216                        B = A(:,p) ;
217                        B = spqr_qmult (N.Q,B,3) ;
218                        B = B*N.X ;
219                    else
220                        % B = A * ( P * ( Q * X ) )
221                        B = spqr_qmult (N.Q, N.X, 1) ;
222                        B = B(N.P,:) ;
223                        B = A * B ;
224                    end
225
226                otherwise
227                    error ('unrecognized N struct') ;
228            end
229
230        otherwise
231            error ('unrecognized method') ;
232    end
233
234else
235
236    %---------------------------------------------------------------------------
237    % N held as an explicit matrix
238    %---------------------------------------------------------------------------
239
240    switch method
241        case 0
242            B = N'*A ;
243        case 1
244            B = N*A ;
245        case 2
246            B = A*N' ;
247        case 3
248            B = A*N ;
249        otherwise
250            error ('unrecognized method') ;
251    end
252
253end
254