Mercurial > octave

comparison libinterp/corefcn/sparse-xpow.cc @ 31235:b542b88ad3b6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use space between function name and '(' in sparse-xpow.cc. * sparse-xpow.cc: Use space between function name and '(' in sparse-xpow.cc.
author Rik <rik@octave.org>
date Tue, 20 Sep 2022 14:43:56 -0700
parents a1318deb4584
children e88a07dec498
comparison
equal deleted inserted replaced
31234:a1318deb4584 31235:b542b88ad3b6
65 65
66 octave_idx_type nr = a.rows (); 66 octave_idx_type nr = a.rows ();
67 octave_idx_type nc = a.cols (); 67 octave_idx_type nc = a.cols ();
68 68
69 if (nr == 0 || nc == 0) 69 if (nr == 0 || nc == 0)
70 return SparseMatrix(); 70 return SparseMatrix ();
71 71
72 // If we are here, A is not empty ==> A needs to be square. 72 // If we are here, A is not empty ==> A needs to be square.
73 if (nr != nc) 73 if (nr != nc)
74 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); 74 error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
75 75
110 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); 110 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
111 } 111 }
112 else 112 else
113 atmp = a; 113 atmp = a;
114 114
115 if (atmp.nnz() == 0) // Fast return for all-zeros matrix 115 if (atmp.nnz () == 0) // Fast return for all-zeros matrix
116 return atmp; 116 return atmp;
117 117
118 SparseMatrix result (atmp); 118 SparseMatrix result (atmp);
119 119
120 btmp--; 120 btmp--;
122 // There are two approaches to the actual exponentiation. 122 // There are two approaches to the actual exponentiation.
123 // Exponentiation by squaring uses only a logarithmic number 123 // Exponentiation by squaring uses only a logarithmic number
124 // of multiplications but the matrices it multiplies tend to be dense 124 // of multiplications but the matrices it multiplies tend to be dense
125 // towards the end. 125 // towards the end.
126 // Linear multiplication uses a linear number of multiplications 126 // Linear multiplication uses a linear number of multiplications
127 // but one of the matrices it uses will be as sparse as the original matrix. 127 // but one of the matrices it uses will be as sparse as the original
128 // matrix.
128 // 129 //
129 // The time to multiply fixed-size matrices is strongly affected by their 130 // The time to multiply fixed-size matrices is strongly affected by their
130 // sparsity. Denser matrices take much longer to multiply together. 131 // sparsity. Denser matrices take much longer to multiply together.
131 // See this URL for a worked-through example: 132 // See this URL for a worked-through example:
132 // https://octave.discourse.group/t/3216/4 133 // https://octave.discourse.group/t/3216/4
133 // 134 //
134 // The tradeoff is between many fast multiplications or a few slow ones. 135 // The tradeoff is between many fast multiplications or a few slow ones.
135 // 136 //
136 // Large exponents favor the squaring technique, and sparse matrices favor 137 // Large exponents favor the squaring technique, and sparse matrices
137 // linear multiplication. 138 // favor linear multiplication.
138 // 139 //
139 // We calculate a threshold based on the sparsity of the input 140 // We calculate a threshold based on the sparsity of the input
140 // and use squaring for exponents larger than that. 141 // and use squaring for exponents larger than that.
141 // 142 //
142 // FIXME: Improve this threshold calculation. 143 // FIXME: Improve this threshold calculation.
143 144
144 uint64_t sparsity = atmp.numel() / atmp.nnz(); // reciprocal of density 145 uint64_t sparsity = atmp.numel () / atmp.nnz (); // reciprocal of density
145 int threshold = (sparsity >= 1000) ? 40 146 int threshold = (sparsity >= 1000) ? 40
146 : (sparsity >= 100) ? 20 147 : (sparsity >= 100) ? 20
147 : 3; 148 : 3;
148 149
149 if (btmp > threshold) // use squaring technique 150 if (btmp > threshold) // use squaring technique
178 179
179 octave_idx_type nr = a.rows (); 180 octave_idx_type nr = a.rows ();
180 octave_idx_type nc = a.cols (); 181 octave_idx_type nc = a.cols ();
181 182
182 if (nr == 0 || nc == 0) 183 if (nr == 0 || nc == 0)
183 return SparseMatrix(); 184 return SparseMatrix ();
184 185
185 // If we are here, A is not empty ==> A needs to be square. 186 // If we are here, A is not empty ==> A needs to be square.
186 if (nr != nc) 187 if (nr != nc)
187 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); 188 error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
188 189
220 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); 221 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
221 } 222 }
222 else 223 else
223 atmp = a; 224 atmp = a;
224 225
225 if (atmp.nnz() == 0) // Fast return for all-zeros matrix 226 if (atmp.nnz () == 0) // Fast return for all-zeros matrix
226 return atmp; 227 return atmp;
227 228
228 SparseComplexMatrix result (atmp); 229 SparseComplexMatrix result (atmp);
229 230
230 btmp--; 231 btmp--;
233 // See the long comment in xpow (const SparseMatrix& a, double b) 234 // See the long comment in xpow (const SparseMatrix& a, double b)
234 // for more details. 235 // for more details.
235 // 236 //
236 // FIXME: Improve this threshold calculation. 237 // FIXME: Improve this threshold calculation.
237 238
238 uint64_t sparsity = atmp.numel() / atmp.nnz(); // reciprocal of density 239 uint64_t sparsity = atmp.numel () / atmp.nnz (); // reciprocal of density
239 int threshold = (sparsity >= 1000) ? 40 240 int threshold = (sparsity >= 1000) ? 40
240 : (sparsity >= 100) ? 20 241 : (sparsity >= 100) ? 20
241 : 3; 242 : 3;
242 243
243 if (btmp > threshold) // use squaring technique 244 if (btmp > threshold) // use squaring technique