1/* linbox/algorithms/lanczos.inl 2 * Copyright (C) 2002 Bradford Hovinen 3 * 4 * Written by Bradford Hovinen <hovinen@cis.udel.edu> 5 * 6 * ------------------------------------ 7 * 8 * 9 * ========LICENCE======== 10 * This file is part of the library LinBox. 11 * 12 * LinBox is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU Lesser General Public 14 * License as published by the Free Software Foundation; either 15 * version 2.1 of the License, or (at your option) any later version. 16 * 17 * This library is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 20 * Lesser General Public License for more details. 21 * 22 * You should have received a copy of the GNU Lesser General Public 23 * License along with this library; if not, write to the Free Software 24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 25 * ========LICENCE======== 26 *. 27 */ 28 29#ifndef __LINBOX_lanczos_INL 30#define __LINBOX_lanczos_INL 31 32#include <vector> 33#include <algorithm> 34 35#include "linbox/blackbox/archetype.h" 36#include "linbox/blackbox/diagonal.h" 37#include "linbox/blackbox/transpose.h" 38#include "linbox/util/debug.h" 39#include "linbox/vector/vector-domain.h" 40 41namespace LinBox 42{ 43 44#ifdef DETAILED_TRACE 45 46 template <class Field, class LVector> 47 void traceReport (std::ostream &out, VectorDomain<Field> &VD, const char *text, size_t iter, const LVector &v) 48 { 49 out << text << " [" << iter << "]: "; 50 VD.write (out, v) << std::endl; 51 } 52 53 template <class Field, class LVector> 54 void traceReport (std::ostream &out, const Field &F, const char *text, size_t iter, const typename Field::Element &a) 55 { 56 out << text << " [" << iter << "]: "; 57 F.write (out, a) << std::endl; 58 } 59 60#else 61 62 template <class Field, class LVector> 63 inline void traceReport (std::ostream &out, VectorDomain<Field> &VD, const char *text, size_t iter, const LVector &v) 64 {} 65 66 template <class Field, class LVector> 67 void traceReport (std::ostream &out, const Field &F, const char *text, size_t iter, const typename Field::Element &a) 68 {} 69 70#endif 71 72 template <class Field, class LVector> 73 template <class Blackbox> 74 LVector &LanczosSolver<Field, LVector>::solve (const Blackbox &A, LVector &x, const LVector &b) 75 { 76 linbox_check ((x.size () == A.coldim ()) && 77 (b.size () == A.rowdim ())); 78 79 commentator().start ("Solving linear system (Lanczos)", "LanczosSolver::solve"); 80 81 bool success = false; 82 LVector d1, d2, b1, b2, bp, y, Ax, ATAx, ATb; 83 84 VectorWrapper::ensureDim (_w[0], A.coldim ()); 85 VectorWrapper::ensureDim (_w[1], A.coldim ()); 86 VectorWrapper::ensureDim (_Aw, A.coldim ()); 87 88 Givaro::GeneralRingNonZeroRandIter<Field> real_ri (_randiter); 89 RandomDenseStream<Field, LVector, Givaro::GeneralRingNonZeroRandIter<Field> > stream (field(), real_ri, A.coldim ()); 90 91 for (unsigned int i = 0; !success && i < _traits.trialsBeforeFailure; ++i) { 92 std::ostream &report = commentator().report (Commentator::LEVEL_UNIMPORTANT, INTERNAL_DESCRIPTION); 93 94 switch (_traits.preconditioner ) { 95 case Preconditioner::None: 96 success = iterate (A, x, b); 97 break; 98 99 case Preconditioner::Symmetrize: 100 { 101 VectorWrapper::ensureDim (bp, A.coldim ()); 102 103 Transpose<Blackbox> AT (&A); 104 Compose<Transpose<Blackbox>, Blackbox> B (&AT, &A); 105 106 AT.apply (bp, b); 107 108 success = iterate (B, x, bp); 109 110 break; 111 } 112 113 case Preconditioner::PartialDiagonal: 114 { 115 VectorWrapper::ensureDim (d1, A.coldim ()); 116 VectorWrapper::ensureDim (y, A.coldim ()); 117 118 stream >> d1; 119 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> D (d1); 120 Compose<Blackbox, Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> > B (&A, &D); 121 122 report << "Random D: "; 123 _VD.write (report, d1) << std::endl; 124 125 if ((success = iterate (B, y, b))) 126 D.apply (x, y); 127 128 break; 129 } 130 131 case Preconditioner::PartialDiagonalSymmetrize: 132 { 133 VectorWrapper::ensureDim (d1, A.rowdim ()); 134 VectorWrapper::ensureDim (b1, A.rowdim ()); 135 VectorWrapper::ensureDim (bp, A.coldim ()); 136 137 stream >> d1; 138 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> D (d1); 139 Transpose<Blackbox> AT (&A); 140 Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, Blackbox> B1 (&D, &A); 141 Compose<Transpose<Blackbox>, Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, Blackbox> > B (&AT, &B1); 142 143 report << "Random D: "; 144 _VD.write (report, d1) << std::endl; 145 146 D.apply (b1, b); 147 AT.apply (bp, b1); 148 149 success = iterate (B, x, bp); 150 151 break; 152 } 153 154 case Preconditioner::FullDiagonal: 155 { 156 VectorWrapper::ensureDim (d1, A.coldim ()); 157 VectorWrapper::ensureDim (d2, A.rowdim ()); 158 VectorWrapper::ensureDim (b1, A.rowdim ()); 159 VectorWrapper::ensureDim (b2, A.coldim ()); 160 VectorWrapper::ensureDim (bp, A.coldim ()); 161 VectorWrapper::ensureDim (y, A.coldim ()); 162 163 stream >> d1 >> d2; 164 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> D1 (d1); 165 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> D2 (d2); 166 Transpose<Blackbox> AT (&A); 167 168 Compose<Blackbox, 169 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> > B1 (&A, &D1); 170 171 Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, 172 Compose<Blackbox, 173 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> > > B2 (&D2, &B1); 174 175 Compose<Transpose<Blackbox>, 176 Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, 177 Compose<Blackbox, 178 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> > > > B3 (&AT, &B2); 179 180 Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, 181 Compose<Transpose<Blackbox>, 182 Compose<Diagonal<Field, typename VectorTraits<LVector>::VectorCategory>, 183 Compose<Blackbox, 184 Diagonal<Field, typename VectorTraits<LVector>::VectorCategory> > > > > B (&D1, &B3); 185 186 report << "Random D_1: "; 187 _VD.write (report, d1) << std::endl; 188 189 report << "Random D_2: "; 190 _VD.write (report, d2) << std::endl; 191 192 D2.apply (b1, b); 193 AT.apply (b2, b1); 194 D1.apply (bp, b2); 195 196 if ((success = iterate (B, y, bp))) 197 D1.apply (x, y); 198 199 break; 200 } 201 202 default: 203 throw PreconditionFailed (__func__, __LINE__, 204 "preconditioner should be None, Symmetrize, PartialDiagonal," 205 "PartialDiagonalSymmetrize, or FullDiagonal"); 206 } 207 208 if (success && _traits.checkResult) { 209 VectorWrapper::ensureDim (Ax, A.rowdim ()); 210 211 if ((_traits.preconditioner == Preconditioner::Symmetrize) || 212 (_traits.preconditioner == Preconditioner::PartialDiagonalSymmetrize) || 213 (_traits.preconditioner == Preconditioner::FullDiagonal)) 214 { 215 VectorWrapper::ensureDim (ATAx, A.coldim ()); 216 VectorWrapper::ensureDim (ATb, A.coldim ()); 217 218 commentator().start ("Checking whether A^T Ax = A^T b"); 219 220 A.apply (Ax, x); 221 A.applyTranspose (ATAx, Ax); 222 A.applyTranspose (ATb, b); 223 224 if (_VD.areEqual (ATAx, ATb)) 225 commentator().stop ("passed"); 226 else { 227 commentator().stop ("FAILED"); 228 success = false; 229 } 230 } 231 else { 232 commentator().start ("Checking whether Ax=b"); 233 234 A.apply (Ax, x); 235 236 if (_VD.areEqual (Ax, b)) 237 commentator().stop ("passed"); 238 else { 239 commentator().stop ("FAILED"); 240 success = false; 241 } 242 } 243 } 244 } 245 246 if (success) { 247 commentator().stop ("done", "Solve successful", "BlockLanczosSolver::solve"); 248 return x; 249 } 250 else { 251 commentator().stop ("done", "Solve failed", "BlockLanczosSolver::solve"); 252 throw SolveFailed (); 253 } 254 } 255 256 template <class Field, class LVector> 257 template<class Blackbox> 258 bool LanczosSolver<Field, LVector>::iterate (const Blackbox &A, LVector &x, const LVector &b) 259 { 260 commentator().start ("Lanczos iteration", "LanczosSolver::iterate", A.coldim ()); 261 262 // j is really a flip-flop: 0 means "even" and 1 means "odd". So "j" and 263 // "j-2" are accessed with [j], while "j-1" and "j+1" are accessed via 264 // [1-j] 265 266 unsigned int j = 1, prods = 1, iter = 2; 267 268 // N.B. For purposes of efficiency, I am purposefully making the 269 // definitions of alpha and beta to be the *negatives* of what are given 270 // in the Lambert thesis. This allows me to use stock vector AXPY 271 // without any special modifications. 272 273 typename Field::Element alpha, beta, delta[2], wb; 274 275 // Zero out the vector _w[0] 276 _VD.subin (_w[0], _w[0]); 277 278 // Get a random vector _w[1] 279 RandomDenseStream<Field, LVector> stream (field(), _randiter, A.coldim ()); 280 stream >> _w[1]; 281 282 std::ostream &report = commentator().report (Commentator::LEVEL_UNIMPORTANT, INTERNAL_DESCRIPTION); 283 284 traceReport (report, _VD, "w", 1, _w[1]); 285 286 A.apply (_Aw, _w[j]); // Aw_j 287 _VD.dot (delta[j], _w[j], _Aw); // delta_j <- <w_j, Aw_j> 288 289 if (field().isZero (delta[j])) { 290 commentator().stop ("FAILED", "<w_1, Aw_1> = 0", "LanczosSolver::iterate"); 291 return false; 292 } 293 294 _VD.dot (alpha, _Aw, _Aw); // alpha <- -<Aw_j, Aw_j> / delta_j 295 field().divin (alpha, delta[j]); 296 field().negin (alpha); 297 298 field().subin (beta, beta); // beta <- 0 299 300 _VD.dot (wb, _w[j], b); // x <- <w_j, b> / delta_j w_j 301 field().divin (wb, delta[j]); 302 _VD.mul (x, _w[j], wb); 303 304 while (!field().isZero (delta[j])) { 305 commentator().progress (); 306 307 report << "Total matrix-vector products so far: " << prods << std::endl; 308 309 // traceReport (report, field(), "alpha", iter, alpha); 310 // traceReport (report, field(), "beta", iter, alpha); 311 traceReport (report, _VD, "w", iter - 1, _w[1 - j]); 312 traceReport (report, _VD, "w", iter, _w[j]); 313 314 _VD.mulin (_w[1 - j], beta); // w_j+1 <- Aw_j + alpha w_j + beta w_j-1 315 _VD.axpyin (_w[1 - j], alpha, _w[j]); 316 _VD.addin (_w[1 - j], _Aw); 317 318 traceReport (report, _VD, "w", iter + 1, _w[1 - j]); 319 traceReport (report, _VD, "Aw", iter, _Aw); 320 321 j = 1 - j; // j <- j + 1 322 323 A.apply (_Aw, _w[j]); // Aw_j 324 325 _VD.dot (delta[j], _w[j], _Aw); // delta_j <- <w_j, Aw_j> 326 327 // traceReport (report, field(), "delta", iter - 1, delta[1 - j]); 328 // traceReport (report, field(), "delta", iter, delta[j]); 329 330 if (!field().isZero (delta[j])) { 331 _VD.dot (alpha, _Aw, _Aw); // alpha <- -<Aw_j, Aw_j> / delta_j 332 field().divin (alpha, delta[j]); 333 field().negin (alpha); 334 335 field().div (beta, delta[j], delta[1 - j]); // beta <- -delta_j / delta_j-1 336 field().negin (beta); 337 338 _VD.dot (wb, _w[j], b); // x <- x + <w_j, b> / delta_j w_j 339 field().divin (wb, delta[j]); 340 _VD.axpyin (x, wb, _w[j]); 341 } 342 343 ++prods; 344 ++iter; 345 } 346 347 commentator().indent (report); 348 report << "Total matrix-vector products: " << prods << std::endl; 349 350 commentator().stop ("done", "delta_j = 0", "LanczosSolver::iterate"); 351 return true; 352 } 353 354} // namespace LinBox 355 356#endif // __LINBOX_lanczos_INL 357 358// Local Variables: 359// mode: C++ 360// tab-width: 4 361// indent-tabs-mode: nil 362// c-basic-offset: 4 363// End: 364// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 365