Ifpack Package Browser (Single Doxygen Collection) Development
Ifpack_SupportGraph.h
Go to the documentation of this file.
00001 /*@HEADER 
00002 // ***********************************************************************
00003 //
00004 //       Ifpack: Object-Oriented Algebraic Preconditioner Package
00005 //                 Copyright (2002) Sandia Corporation
00006 //                                        
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //                                                  
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 // 
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of       
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU   
00018 // Lesser General Public License for more details.     
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //     
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #ifndef IFPACK_SUPPORTGRAPH_H
00031 #define IFPACK_SUPPORTGRAPH_H
00032 
00033 #include "Ifpack_ConfigDefs.h"
00034 #include "Ifpack_Condest.h"
00035 #include "Ifpack_Preconditioner.h"
00036 #include "Ifpack_Amesos.h"
00037 #include "Ifpack_Condest.h"
00038 #include "Epetra_MultiVector.h"
00039 #include "Epetra_Map.h"
00040 #include "Epetra_Comm.h"
00041 #include "Epetra_Time.h"
00042 #include "Epetra_LinearProblem.h"
00043 #include "Epetra_RowMatrix.h"
00044 #include "Epetra_CrsMatrix.h"
00045 
00046 #include "Teuchos_ParameterList.hpp"
00047 #include "Teuchos_RefCountPtr.hpp"
00048 
00049 #include <boost/graph/adjacency_list.hpp>
00050 #include <boost/graph/kruskal_min_spanning_tree.hpp>
00051 #include <boost/graph/prim_minimum_spanning_tree.hpp>
00052 #include <boost/config.hpp>
00053 
00054 using Teuchos::RefCountPtr;
00055 using Teuchos::rcp;
00056 typedef std::pair<int, int> E;
00057 using namespace boost;
00058 
00059 typedef adjacency_list < vecS, vecS, undirectedS,
00060   no_property, property < edge_weight_t, double > > Graph;
00061 typedef graph_traits < Graph >::edge_descriptor Edge;
00062 typedef graph_traits < Graph >::vertex_descriptor Vertex;
00063 
00064 
00065 
00066 template<typename T=Ifpack_Amesos> class Ifpack_SupportGraph : 
00067 public virtual Ifpack_Preconditioner 
00068 {
00069 
00070  public:
00071  
00073 
00075  Ifpack_SupportGraph(Epetra_RowMatrix* Matrix_in);
00076 
00078 
00079 
00081 
00082 
00091  virtual int SetUseTranspose(bool UseTranspose_in);
00092 
00094  
00095 
00097  
00099 
00107  virtual int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00108 
00110 
00121  virtual int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00122 
00124  virtual double NormInf() const {return(0.0);}
00125 
00127 
00128 
00130 
00132  virtual const char * Label() const;
00133 
00135  virtual bool UseTranspose() const {return(UseTranspose_);};
00136 
00138  virtual bool HasNormInf() const {return(false);};
00139 
00141  virtual const Epetra_Comm & Comm() const {return(Matrix_->Comm());};
00142 
00144  virtual const Epetra_Map & OperatorDomainMap() const {return(Matrix_->OperatorDomainMap());};
00145 
00147  virtual const Epetra_Map & OperatorRangeMap() const {return(Matrix_->OperatorRangeMap());};
00148 
00150 
00151 
00153 
00155  virtual bool IsInitialized() const
00156  {
00157    return(IsInitialized_);
00158  }
00159 
00161  virtual bool IsComputed() const
00162  {
00163    return(IsComputed_);
00164  }
00165 
00167 
00179  virtual int SetParameters(Teuchos::ParameterList& List);
00180 
00182 
00185  virtual int Initialize();
00187 
00190  virtual int Compute();
00191 
00193 
00194 
00196 
00197 
00199 
00202  virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
00203       const int MaxIters = 1550,
00204       const double Tol = 1e-9,
00205       Epetra_RowMatrix* Matrix_in = 0);
00206 
00208  virtual double Condest() const
00209  {
00210    return(Condest_);
00211  }
00212 
00214  virtual const Epetra_RowMatrix& Matrix() const
00215  {
00216    return(*Matrix_);
00217  }
00218 
00220  virtual std::ostream& Print(std::ostream&) const;
00221 
00223  virtual int NumInitialize() const
00224  {
00225    return(NumInitialize_);
00226  }
00227 
00229  virtual int NumCompute() const
00230  {
00231    return(NumCompute_);
00232  }
00233 
00235  virtual int NumApplyInverse() const
00236  {
00237    return(NumApplyInverse_);
00238  }
00239 
00241  virtual double InitializeTime() const
00242  {
00243    return(InitializeTime_);
00244  }
00245 
00247  virtual double ComputeTime() const
00248  {
00249    return(ComputeTime_);
00250  }
00251 
00253  virtual double ApplyInverseTime() const
00254  {
00255    return(ApplyInverseTime_);
00256  }
00257 
00259  virtual double InitializeFlops() const
00260  {
00261    return(InitializeFlops_);
00262  }
00263 
00265  virtual double ComputeFlops() const
00266  {
00267    return(ComputeFlops_);
00268  }
00269 
00271  virtual double ApplyInverseFlops() const
00272  {
00273    return(ApplyInverseFlops_);
00274  }
00275 
00276 
00278 
00279  protected:
00280 
00281  // Finds the size of a subtree rooted at a vertex (work in progress)
00282  int treecount(const std::vector<Vertex>& v, int *subtreesize, int node);
00283 
00284  // Partitions the spanning tree (work in progress)
00285  int treepartition(int *table, int* children, int *subtreesize, std::vector<int>& roots, 
00286        int node, int n, int t);
00287 
00288  // Finds the largest edge between two trees (work in progress)
00289  int largestbetween(int *table, int* children, const Graph& graph, const property_map<Graph, edge_weight_t>::type& map,
00290         int tree1, int tree2, double *largest, int *extrasource, int *extratarget, int num_verts);
00291 
00292  // Finds a list of all the vertices of a tree rooted at input vertex (work in progress)
00293  int findall(std::vector<int>& tree, int root, int *table, int *children, int num_verts);
00294 
00295 
00297  int AMST();
00298  
00300  int FindSupport();
00301 
00303  Teuchos::RefCountPtr<const Epetra_RowMatrix> Matrix_;
00304  
00306  Teuchos::RefCountPtr<Epetra_CrsMatrix> Support_;
00307 
00309  string Label_;
00310 
00312  bool IsInitialized_;
00313 
00315  bool IsComputed_;
00316 
00318  bool UseTranspose_;
00319 
00321  Teuchos::ParameterList List_;
00322 
00324  double Condest_;
00325 
00327  int NumInitialize_;
00328  
00330  int NumCompute_;
00331  
00333  mutable int NumApplyInverse_;
00334  
00336  double InitializeTime_;
00337 
00339  double ComputeTime_;
00340  
00342  mutable double ApplyInverseTime_;
00343 
00345  double InitializeFlops_;
00346 
00348  double ComputeFlops_;
00349 
00351  mutable double ApplyInverseFlops_;
00352  
00354  Teuchos::RefCountPtr<Epetra_Time> Time_;
00355 
00357  Teuchos::RefCountPtr<T> Inverse_;
00358 
00360  int NumForests_;
00361 
00363  double Offset_;
00364 
00366  int KeepDiag_;
00367 
00368 }; // class Ifpack_SupportGraph<T>
00369 
00370 
00371 
00372 //==============================================================================
00373 template<typename T>
00374 Ifpack_SupportGraph<T>::Ifpack_SupportGraph(Epetra_RowMatrix* Matrix_in):
00375 Matrix_(rcp(Matrix_in,false)),
00376   IsInitialized_(false),
00377   IsComputed_(false),
00378   UseTranspose_(false),
00379   Condest_(-1.0),
00380   NumInitialize_(0),
00381   NumCompute_(0),
00382   NumApplyInverse_(0),
00383   InitializeTime_(0.0),
00384   ComputeTime_(0.0),
00385   ApplyInverseTime_(0.0),
00386   InitializeFlops_(0.0),
00387   ComputeFlops_(0.0),
00388   ApplyInverseFlops_(0.0),
00389   NumForests_(1),
00390   Offset_(1),
00391   KeepDiag_(1)
00392 {
00393   
00394   Teuchos::ParameterList List_in;
00395   SetParameters(List_in);
00396 }
00397 //============================================================================== 
00398 template<typename T>
00399 int Ifpack_SupportGraph<T>::treecount(const std::vector<Vertex>& v, int *subtreesize, int node)
00400 {
00401   int parent = v[node];
00402  
00403   if(node != parent)
00404     {
00405       subtreesize[node] = subtreesize[node] + 1;
00406       treecount(v, subtreesize, parent);
00407     }
00408 
00409   return 0;
00410 }
00411 //==============================================================================
00412 template<typename T>
00413 int Ifpack_SupportGraph<T>::treepartition(int *table, int* children, int *subtreesize, 
00414             std::vector<int>& roots,int node, int n, int t)
00415 {
00416   subtreesize[node] = 1;
00417   for(int i = table[node]; i < table[node+1]; i++)
00418     {
00419       int child = children[i];
00420     
00421       if(subtreesize[child] > (n/t + 1))
00422   {
00423     treepartition(table, children, subtreesize, roots, child, n, t);
00424   }
00425 
00426       if(subtreesize[child] > (n/t))
00427   {
00428     children[i] = -children[i];
00429     roots.push_back(child);
00430     
00431   }
00432       else
00433   {
00434     subtreesize[node] = subtreesize[node] + subtreesize[child];
00435   }
00436     }
00437   return 0;
00438 }
00439 //==============================================================================
00440 template<typename T>
00441 int Ifpack_SupportGraph<T>::largestbetween(int* table, int* children, const Graph& graph, 
00442               const property_map<Graph, edge_weight_t>::type& map, 
00443               int tree1, int tree2,
00444              double *largest, int *extrasource, int *extratarget, int num_verts)
00445 {
00446   std::vector <int> subtree1;
00447   std::vector <int> subtree2;
00448 
00449   if(tree1 < 0)
00450     {
00451       tree1 = -tree1;
00452     }
00453 
00454   if(tree2 < 0)
00455     {
00456       tree2 = -tree2;
00457     }
00458 
00459   subtree1.push_back(tree1);
00460   subtree2.push_back(tree2);
00461 
00462 
00463   findall(subtree1, tree1, table, children, num_verts);
00464   findall(subtree2, tree2, table, children, num_verts);
00465 
00466 
00467   for(int i = 0; i < subtree1.size(); i++)
00468     {
00469       for(int j = 0; j < subtree2.size(); j++)
00470   {
00471 
00472     if(edge(subtree1[i], subtree2[j], graph).second)
00473       {                                                  
00474 
00475         double temp = get(map, edge(subtree1[i],subtree2[j], graph).first);                                                                                  
00476 
00477         if(temp < *largest)         
00478     {
00479       *largest = temp;
00480       *extrasource = subtree1[i];
00481       *extratarget = subtree2[j];
00482     }
00483 
00484       }    
00485   }
00486     }
00487 
00488 
00489 
00490     return 0;
00491 }
00492 //==============================================================================
00493 template<typename T>
00494 int Ifpack_SupportGraph<T>::findall(std::vector<int>& tree, int root, int *table, int *children, int num_verts)
00495 {
00496   int upper;
00497   if(root == num_verts-1)
00498     {
00499       upper = num_verts;
00500     }
00501   else
00502     {
00503       upper = table[root+1];
00504     }
00505   for(int i = table[root]; i < upper; i++)
00506     {
00507       int child = children[i];
00508      
00509       if(child > 0)
00510   {
00511     tree.push_back(child);
00512     findall(tree, child, table, children,num_verts);
00513   }
00514     }
00515 }
00516 //============================================================================== 
00517 template<typename T>
00518 int Ifpack_SupportGraph<T>::AMST()
00519 {
00520   
00521  
00522 
00523   // Extract matrix dimensions
00524   int rows = (*Matrix_).NumGlobalRows();
00525   int cols = (*Matrix_).NumGlobalCols();
00526   int num_edges  = ((*Matrix_).NumMyNonzeros() - (*Matrix_).NumMyDiagonals())/2;
00527 
00528   // Assert square matrix
00529   IFPACK_CHK_ERR((rows == cols));
00530 
00531   // Rename for clarity
00532   int num_verts = rows;
00533 
00534   double fill = .9;
00535   int t = fill*pow(num_verts,.5);
00536 
00537   // Create data structures for the BGL code and temp data structures for extraction 
00538   E *edge_array = new E[num_edges];
00539   double *weights = new double[num_edges];
00540   double *shiftedweights = new double[num_edges];
00541 
00542   int num_entries;
00543   int max_num_entries = (*Matrix_).MaxNumEntries();
00544   double *values = new double[max_num_entries];
00545   int *indices = new int[max_num_entries];
00546 
00547   double * diagonal = new double[num_verts];
00548   double shift = 0;
00549 
00550   for(int i = 0; i < max_num_entries; i++)
00551     {
00552       values[i]=0;
00553       indices[i]=0;
00554     }
00555 
00556   // Extract from the epetra matrix keeping only one edge per pair (assume symmetric) 
00557   int k = 0;
00558   for(int i = 0; i < num_verts; i++)
00559     {
00560       (*Matrix_).ExtractMyRowCopy(i,max_num_entries,num_entries,values,indices);
00561 
00562       for(int j = 0; j < num_entries; j++)
00563         {
00564 
00565           if(i == indices[j])
00566             {
00567               diagonal[i] = values[j];
00568             }
00569 
00570           if(i < indices[j])
00571             {
00572               edge_array[k] = E(i,indices[j]);
00573               if(values[j] < shift)
00574                 shift = values[j];
00575 
00576               weights[k] = values[j];
00577 
00578               k++;
00579             }
00580         }
00581     }
00582 
00583   shift = shift - 1;
00584 
00585   for(int i = 0; i < num_edges; i++)
00586     {
00587       shiftedweights[i] = weights[i] - shift;
00588     }
00589 
00590 
00591   std::vector<int> TreeNz(num_verts,1);
00592 
00593   int *TotalNz = new int[num_verts];
00594   for(int i = 0; i < num_verts; i++)
00595     {
00596       TotalNz[i] = 1;
00597     }
00598 
00599 
00600   Graph gtemp(edge_array, edge_array + num_edges, shiftedweights, num_verts);
00601   //gtemp = Graph(edge_array, edge_array + num_edges, shiftedweights, num_verts);
00602 
00603   property_map<Graph, edge_weight_t>::type weightmap = get(edge_weight, gtemp);
00604 
00605   //weightmap = get(edge_weight, gtemp);
00606 
00607   std::vector < Vertex > p(num_vertices(gtemp));
00608   prim_minimum_spanning_tree(gtemp, &p[0]);
00609 
00610   std::vector<int> roots;
00611   int numchildren[num_verts];
00612   int table[num_verts];
00613   int children[num_verts];
00614   int *subtreesize = new int[num_verts];
00615   for(std::size_t i = 0; i != p.size(); ++i)
00616     {
00617       numchildren[i] = 0;
00618       table[i] = 0;
00619       children[i] = 0;
00620       subtreesize[i] = 0;
00621     }
00622 
00623   for (std::size_t i = 0; i != p.size(); ++i)
00624     {
00625 
00626       if (p[i] != i)
00627   {
00628     numchildren[p[i]] = numchildren[p[i]] + 1;
00629     TreeNz[p[i]] = TreeNz[p[i]] + 1;
00630     TreeNz[i] = TreeNz[i] + 1;
00631     TotalNz[p[i]] = TotalNz[p[i]] + 1;
00632     TotalNz[i] = TotalNz[i] + 1;
00633     //std::cout << "parent[" << i << "] = " << p[i] << std::endl;
00634   }
00635       else
00636   {
00637     roots.push_back(i);
00638   }
00639     }
00640 
00641 
00642 
00643 
00644   k = 0;
00645  
00646   for (std::size_t i = 0; i != p.size(); ++i)
00647     {
00648       table[i] = k;
00649       k = k + numchildren[i];
00650     }
00651  
00652  
00653   for (std::size_t i = 0; i != p.size(); ++i)
00654     {
00655       if(i != p[i])
00656   {
00657     //    std::cout << table[p[i]] << "    " << numchildren[p[i]] << "     " << table[p[i]] + numchildren[i] - 1 << std::endl;
00658     children[table[p[i]] + numchildren[p[i]] - 1] = i;
00659     numchildren[p[i]] = numchildren[p[i]] - 1;
00660   }
00661     }
00662 
00663  
00664   for(std::size_t i = 0; i != p.size(); ++i)
00665     {
00666       treecount(p,subtreesize, i);
00667     }
00668 
00669  
00670   for(int i = 0; i < roots.size(); i++)
00671     {
00672       treepartition(table, children, subtreesize, roots, roots[i], num_verts, t);
00673     }
00674 
00675   /*
00676   for(int i = 0; i < roots.size(); i++)
00677     {
00678       p[roots[i]] = i;
00679       std::cout << roots[i] << std::endl;
00680       }*/
00681   /*
00682   std::cout << "children" << std::endl;
00683   for(std::size_t i = 0; i != p.size(); ++i)
00684     {
00685       std::cout << children[i] << std::endl;
00686       }*/
00687 
00688   std::vector<int> ExtraIndices[num_verts];
00689   std::vector<double> ExtraValues[num_verts];
00690 
00691   for(int i = 0; i < num_verts; i++)
00692     {
00693       std::vector<int> temp;
00694       std::vector<double> temp2;
00695       ExtraIndices[i] = temp;
00696       ExtraValues[i] = temp2;
00697     }
00698 
00699 
00700 
00701   for(int i = 0; i < roots.size(); i++)
00702     {
00703       for(int j = i+1; j < roots.size(); j++)
00704   {
00705     double largest = -shift;
00706     int extrasource = -1;
00707     int extratarget = -1;
00708 
00709 
00710     largestbetween(table, children, gtemp, weightmap, roots[i],roots[j],&largest,&extrasource,&extratarget,num_verts);
00711 
00712     if(largest < -shift)
00713       {
00714 
00715         if((p[extrasource] != extratarget) && (p[extratarget] != extrasource))
00716     {
00717 
00718       ExtraIndices[extrasource].push_back(extratarget);
00719       ExtraIndices[extratarget].push_back(extrasource);
00720       ExtraValues[extrasource].push_back(largest+shift);
00721       ExtraValues[extratarget].push_back(largest+shift);
00722       TotalNz[extrasource] = TotalNz[extrasource] + 1;
00723       TotalNz[extratarget] = TotalNz[extratarget] + 1;
00724     }
00725       } 
00726     
00727   }
00728     
00729     
00730 
00731     }
00732 
00733 
00734   
00735   Support_ = rcp(new Epetra_CrsMatrix(Copy, Matrix().RowMatrixRowMap(),TotalNz, true));
00736 
00737 
00738   for(int i = 0; i < num_verts; i++)
00739     {
00740 
00741       std::vector<int> Indices(TotalNz[i],0);
00742       std::vector<double> Values(TotalNz[i],0);
00743      
00744       int other;
00745       int upper;
00746       if(p[i] == i)
00747   {
00748     upper = TreeNz[i] - 1;
00749         
00750   }
00751       else
00752   {
00753 
00754     std::cout << TreeNz[i] << std::endl;
00755     upper = TreeNz[i] - 2;
00756     Indices[TreeNz[i]-1] = p[i];
00757 
00758     if(!edge(i,p[i],gtemp).second)
00759       {
00760         std::cout << "WTFFFFFF" << std::endl;
00761       }
00762     Values[TreeNz[i]-1] = get(weightmap, edge(i, p[i], gtemp).first) + shift;
00763   }
00764 
00765       for(int j = 0; j < upper; j++)
00766   {
00767 
00768     other = children[table[i]+j];
00769     if(other < 0)
00770       other = -other;
00771 
00772     Indices[j+1] = other;
00773     Values[j+1] = get(weightmap, edge(i, other, gtemp).first) + shift;
00774   }
00775 
00776       int s = 0;
00777       for(int j = TreeNz[i] + 1; j < TotalNz[i]; j++)
00778   {
00779         std::cout << "extra" << std::endl;
00780         Indices[j] = ExtraIndices[i][s];
00781         Values[j] = ExtraValues[i][s];
00782 
00783         s++;
00784 
00785   }
00786     
00787 
00788       Indices[0] = i;
00789       Values[0] = diagonal[i];
00790 
00791       (*Support_).InsertGlobalValues(i,TotalNz[i],&Values[0],&Indices[0]);
00792 
00793     }
00794 
00795   (*Support_).FillComplete();
00796 
00797   //(*Support_).Print(std::cout);
00798 
00799   delete subtreesize;
00800   delete TotalNz;
00801 
00802   return 0;
00803 }
00804 //============================================================================== 
00805 template<typename T>
00806 int Ifpack_SupportGraph<T>::FindSupport()
00807 {
00808 
00809   // Extract matrix dimensions                                                                  
00810   int rows = (*Matrix_).NumGlobalRows();
00811   int cols = (*Matrix_).NumGlobalCols();
00812   int num_edges  = ((*Matrix_).NumMyNonzeros() - (*Matrix_).NumMyDiagonals())/2;
00813 
00814   // Assert square matrix                                                                       
00815   IFPACK_CHK_ERR((rows == cols));
00816 
00817   // Rename for clarity                                                                         
00818   int num_verts = rows;
00819 
00820   // Create data structures for the BGL code and temp data structures for extraction            
00821   E *edge_array = new E[num_edges];
00822   double *weights = new double[num_edges];
00823  
00824   int num_entries;
00825   int max_num_entries = (*Matrix_).MaxNumEntries();
00826   double *values = new double[max_num_entries];
00827   int *indices = new int[max_num_entries];
00828 
00829   double * diagonal = new double[num_verts];
00830 
00831   
00832   for(int i = 0; i < max_num_entries; i++)
00833     {
00834       values[i]=0;
00835       indices[i]=0;
00836     }
00837 
00838   // Extract from the epetra matrix keeping only one edge per pair (assume symmetric)           
00839   int k = 0;
00840   for(int i = 0; i < num_verts; i++)
00841     {
00842       (*Matrix_).ExtractMyRowCopy(i,max_num_entries,num_entries,values,indices);
00843 
00844       for(int j = 0; j < num_entries; j++)
00845   {
00846 
00847     if(i == indices[j])
00848       {
00849         diagonal[i] = values[j];
00850       }
00851 
00852     if(i < indices[j])
00853       {
00854         edge_array[k] = E(i,indices[j]);
00855         weights[k] = values[j];
00856 
00857         k++;
00858       }
00859   }
00860     }
00861   
00862   // Create BGL graph                                                                           
00863   Graph g(edge_array, edge_array + num_edges, weights, num_verts);
00864   
00865 
00866   property_map < Graph, edge_weight_t >::type weight = get(edge_weight, g);
00867 
00868   std::vector < Edge > spanning_tree;
00869 
00870   // Run Kruskal, actually maximal weight ST since edges are negative                           
00871   kruskal_minimum_spanning_tree(g, std::back_inserter(spanning_tree));
00872   
00873 
00874   std::vector<int> NumNz(num_verts,1);
00875 
00876   //Find the degree of all the vertices                                                         
00877   for (std::vector < Edge >::iterator ei = spanning_tree.begin();
00878        ei != spanning_tree.end(); ++ei)
00879     {
00880       NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
00881       NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
00882     }
00883   
00884   
00885   // Create an stl vector of stl vectors to hold indices and values (neighbour edges)
00886   std::vector< std::vector< int > > Indices(num_verts);
00887   //std::vector<int> Indices[num_verts];
00888   //std::vector<double> Values[num_verts];
00889 
00890   std::vector< std::vector< double > > Values(num_verts);
00891   
00892   for(int i = 0; i < num_verts; i++)
00893     {
00894       std::vector<int> temp(NumNz[i],0);
00895       std::vector<double> temp2(NumNz[i],0);
00896       Indices[i] = temp;
00897       Values[i] = temp2;
00898     }
00899   
00900   int *l = new int[num_verts];
00901   for(int i = 0; i < num_verts; i++)
00902     {
00903       l[i] = 1;
00904     }
00905   
00906   for(int i = 0; i < NumForests_; i++)
00907     {
00908       if(i > 0)
00909   {
00910     spanning_tree.clear();
00911     kruskal_minimum_spanning_tree(g,std::back_inserter(spanning_tree));
00912     for(std::vector < Edge >::iterator ei = spanning_tree.begin();
00913         ei != spanning_tree.end(); ++ei)
00914       {
00915         NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
00916         NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
00917       }
00918     for(int i = 0; i < num_verts; i++)
00919       {
00920         Indices[i].resize(NumNz[i]);
00921         Values[i].resize(NumNz[i]);
00922       }
00923   }
00924 
00925       for (std::vector < Edge >::iterator ei = spanning_tree.begin();
00926      ei != spanning_tree.end(); ++ei)
00927   {
00928     if(KeepDiag_ == 0)
00929       {
00930         Values[source(*ei,g)][0] = Values[source(*ei,g)][0] - weight[*ei];
00931         Values[target(*ei,g)][0] = Values[target(*ei,g)][0] - weight[*ei];
00932       }
00933     Indices[source(*ei,g)][0] = source(*ei,g);
00934     Indices[target(*ei,g)][0] = target(*ei,g);
00935 
00936 
00937     Indices[source(*ei,g)][l[source(*ei,g)]] = target(*ei,g);
00938     Values[source(*ei,g)][l[source(*ei,g)]] = weight[*ei];
00939     l[source(*ei,g)] = l[source(*ei,g)] + 1;
00940 
00941     Indices[target(*ei,g)][l[target(*ei,g)]] = source(*ei,g);
00942     Values[target(*ei,g)][l[target(*ei,g)]] = weight[*ei];
00943     l[target(*ei,g)] = l[target(*ei,g)] + 1;
00944 
00945     remove_edge(*ei,g);
00946   }
00947 
00948     }
00949 
00950   
00951   if(KeepDiag_ == 1)
00952     {
00953       for(int i = 0; i < num_verts; i++)
00954   {
00955     Values[i][0] = diagonal[i];
00956   }
00957     }
00958   
00959   // Create the CrsMatrix for the support graph                                                 
00960   Support_ = rcp(new Epetra_CrsMatrix(Copy, Matrix().RowMatrixRowMap(),l, true));
00961 
00962  
00963   // Fill in the matrix with the stl vectors for each row                                       
00964   for(int i = 0; i < num_verts; i++)
00965     {
00966       Values[i][0] = Values[i][0] + Offset_;
00967 
00968       (*Support_).InsertGlobalValues(i,l[i],&Values[i][0],&Indices[i][0]);
00969     }
00970  
00971   (*Support_).FillComplete();
00972 
00973   //(*Support_).Print(std::cout);    
00974 
00975   delete edge_array;
00976   delete weights;
00977   delete values;
00978   delete indices;
00979   delete l;
00980   delete diagonal;
00981 
00982 
00983   return(0);
00984 }
00985 //==============================================================================
00986 template<typename T>
00987 int Ifpack_SupportGraph<T>::SetParameters(Teuchos::ParameterList& List_in)
00988 {
00989   List_ = List_in;
00990   NumForests_ = List_in.get("MST: forest number", NumForests_);
00991   Offset_ = List_in.get("MST: diagonal offset", Offset_);
00992   KeepDiag_ = List_in.get("MST: keep diagonal", KeepDiag_);
00993 
00994   return(0);
00995 }
00996 //==============================================================================    
00997 template<typename T>
00998 int Ifpack_SupportGraph<T>::Initialize()
00999 {
01000   IsInitialized_ = false;
01001   IsComputed_ = false;
01002 
01003   
01004   
01005   if (Time_ == Teuchos::null)
01006     {
01007       Time_ = Teuchos::rcp( new Epetra_Time(Comm()) );
01008     }
01009 
01010   
01011   Time_->ResetStartTime();
01012  
01013   FindSupport();
01014   //AMST();
01015   Inverse_ = Teuchos::rcp(new T(Support_.get()));
01016 
01017   IFPACK_CHK_ERR(Inverse_->Initialize());
01018 
01019   IsInitialized_ = true;
01020   ++NumInitialize_;
01021   InitializeTime_ += Time_->ElapsedTime();
01022 
01023   return(0);
01024 
01025 }
01026 //==============================================================================
01027 template<typename T>
01028 int Ifpack_SupportGraph<T>::Compute()
01029 {
01030   if (IsInitialized() == false)
01031     IFPACK_CHK_ERR(Initialize());
01032 
01033   Time_->ResetStartTime();
01034   IsComputed_ = false;
01035   Condest_ = -1.0;
01036 
01037   IFPACK_CHK_ERR(Inverse_->Compute());
01038 
01039   IsComputed_ = true;                  
01040   ++NumCompute_;
01041   ComputeTime_ += Time_->ElapsedTime();
01042 
01043 
01044   return(0);
01045 }
01046 //============================================================================== 
01047 template<typename T>
01048 int Ifpack_SupportGraph<T>::SetUseTranspose(bool UseTranspose_in)
01049 {
01050   // store the flag -- it will be set in Initialize() if Inverse_ does not         
01051   // exist.   
01052   UseTranspose_ = UseTranspose_in;
01053 
01054   // If Inverse_ exists, pass it right now.                  
01055   if (Inverse_!=Teuchos::null)
01056     IFPACK_CHK_ERR(Inverse_->SetUseTranspose(UseTranspose_in));
01057   
01058   return(0);
01059 }
01060 //==============================================================================               
01061 template<typename T>
01062 int Ifpack_SupportGraph<T>::
01063 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01064 {
01065   IFPACK_CHK_ERR(Matrix_->Apply(X,Y));
01066   return(0);
01067 }
01068 //==============================================================================                  
01069 template<typename T>
01070 const char * Ifpack_SupportGraph<T>::Label() const
01071 {
01072   return(Label_.c_str());
01073 }
01074 //==============================================================================                  
01075 template<typename T>
01076 int Ifpack_SupportGraph<T>::
01077 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
01078 {
01079   if (!IsComputed())
01080     IFPACK_CHK_ERR(-3);
01081 
01082   Time_->ResetStartTime();
01083 
01084   Inverse_->ApplyInverse(X,Y);
01085 
01086   ++NumApplyInverse_;
01087   ApplyInverseTime_ += Time_->ElapsedTime();
01088 
01089   return(0);
01090 }
01091 //==============================================================================                  
01092 template<typename T>
01093 std::ostream& Ifpack_SupportGraph<T>::
01094 Print(std::ostream& os) const
01095 {
01096   os << "================================================================================" << std::endl;
01097    os << "Ifpack_SupportGraph: " << Label () << endl << endl;
01098   os << "Condition number estimate = " << Condest() << endl;
01099   os << "Global number of rows            = " << Matrix_->NumGlobalRows() << endl;
01100   os << "Number of off diagonal entries in support graph matrix     = " << Support_->NumGlobalNonzeros()-Support_->NumGlobalDiagonals() << endl;
01101   os << "Fraction of off diagonals of support graph/off diagonals of original     = "
01102      << ((double)Support_->NumGlobalNonzeros()-Support_->NumGlobalDiagonals())/(Matrix_->NumGlobalNonzeros()-Matrix_->NumGlobalDiagonals());
01103   os << endl;
01104   os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
01105   os << "-----           -------   --------------       ------------     --------" << endl;
01106   os << "Initialize()    "   << std::setw(10) << NumInitialize_
01107      << "  " << std::setw(15) << InitializeTime_
01108      << "        0.0              0.0" << endl;
01109   os << "Compute()       "   << std::setw(10) << NumCompute_
01110      << "  " << std::setw(22) << ComputeTime_
01111      << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_;
01112   if (ComputeTime_ != 0.0)
01113     os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << endl;
01114   else
01115     os << "     " << std::setw(15) << 0.0 << endl;
01116   os << "ApplyInverse()  "   << std::setw(10) << NumApplyInverse_
01117      << "  " << std::setw(22) << ApplyInverseTime_
01118      << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
01119   if (ApplyInverseTime_ != 0.0)
01120     os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << endl;
01121   else
01122     os << "  " << std::setw(15) << 0.0 << endl;
01123 
01124   os << std::endl << std::endl;
01125   os << "Now calling the underlying preconditioner's print()" << std::endl;
01126 
01127   Inverse_->Print(os);
01128 }
01129 //==============================================================================
01130 template<typename T>
01131 double Ifpack_SupportGraph<T>::
01132 Condest(const Ifpack_CondestType CT, const int MaxIters,
01133   const double Tol, Epetra_RowMatrix* Matrix_in)
01134 {
01135   if (!IsComputed()) // cannot compute right now
01136     {                
01137       return(-1.0);
01138     }
01139  
01140   Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix_in);
01141   
01142   return(Condest_);
01143 }
01144 
01145 #endif // IFPACK_SUPPORTGRAPH_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines