This file is indexed.

/usr/include/trilinos/ml_MultiLevelPreconditioner.h is in libtrilinos-ml-dev 12.10.1-3.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
/*!
 * \file ml_MultiLevelPreconditioner.h
 *
 * \class MultiLevelPreconditioner
 *
 * \brief ML black-box preconditioner for Epetra_RowMatrix derived classes.
 *
 */
/* ******************************************************************** */
/* See the file COPYRIGHT for a complete copyright notice, contact      */
/* person and disclaimer.                                               */
/* ******************************************************************** */
/*#############################################################################
# CVS File Information
#    Current revision: $Revision$
#    Branch:           $Branch$
#    Last modified:    $Date$
#    Modified by:      $Author$
#############################################################################*/

#ifndef ML_MULTILEVELPRECONDITIONER_H
#define ML_MULTILEVELPRECONDITIONER_H

#include "ml_include.h"

#if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS)
// define the following to allow compilation without AztecOO
#ifndef HAVE_ML_AZTECOO
#ifndef AZ_PROC_SIZE
#define AZ_PROC_SIZE 1
#endif
#ifndef AZ_OPTIONS_SIZE
#define AZ_OPTIONS_SIZE 1
#endif
#ifndef AZ_PARAMS_SIZE
#define AZ_PARAMS_SIZE 1
#endif
#ifndef AZ_STATUS_SIZE
#define AZ_STATUS_SIZE 1
#endif
#endif

class Epetra_Map;
class Epetra_BlockMap;
class Epetra_MultiVector;
class Epetra_Comm;
class Epetra_CrsMatrix;
class Epetra_FECrsMatrix;
class Epetra_VbrMatrix;

#include "Epetra_SerialDenseMatrix.h"
#include "Epetra_SerialDenseVector.h"
#include "Epetra_SerialDenseSolver.h"

#define ML_MEM_SIZE      20
#define ML_MEM_INITIAL    0
#define ML_MEM_FINAL      1
#define ML_MEM_SMOOTHER   2
#define ML_MEM_COARSE     3
#define ML_MEM_HIERARCHY  4
#define ML_MEM_PREC_FIRST 5
#define ML_MEM_PREC_OTHER 6
#define ML_MEM_TOT1       7
#define ML_MEM_TOT2       8
#define ML_MEM_INITIAL_MALLOC    10
#define ML_MEM_FINAL_MALLOC      11
#define ML_MEM_SMOOTHER_MALLOC   12
#define ML_MEM_COARSE_MALLOC     13
#define ML_MEM_HIERARCHY_MALLOC  14
#define ML_MEM_PREC_FIRST_MALLOC 15
#define ML_MEM_PREC_OTHER_MALLOC 16
#define ML_MEM_TOT1_MALLOC       17
#define ML_MEM_TOT2_MALLOC       18

#include "Epetra_Operator.h"
#include "Epetra_RowMatrix.h"
#ifdef HAVE_ML_AZTECOO
#include "Epetra_MultiVector.h"
#include "Epetra_MsrMatrix.h"
#endif
#include "Teuchos_ParameterList.hpp"

#ifdef HAVE_ML_EPETRAEXT
#include "EpetraExt_SolverMap_CrsMatrix.h"
#endif
#include "ml_epetra_utils.h"

namespace ML_Epetra
{

  //! Sets default parameters for aggregation-based preconditioners.
  /*! This function, defined in the namespace ML_Epetra, can be used to set
    default values in a user's defined Teuchos::ParameterList.
    \param ProblemType (In) : a std::string, whose possible values are:
       - "SA" : classical smoothed aggregation preconditioners;
       - "NSSA" : default values for Petrov-Galerkin preconditioner for nonsymmetric systems
       - "maxwell" : default values for aggregation preconditioner for eddy current systems
       - "DD" : defaults for 2-level domain decomposition preconditioners based
       on aggregation;
       - "DD-LU" : Like "DD", but use exact LU decompositions on each subdomain;
       - "DD-ML" : 3-level domain decomposition preconditioners, with coarser
       spaces defined by aggregation;
      - "DD-ML-LU" : Like "DD-ML", but with LU decompositions on each subdomain.
    \param List (Out) : list which will populated by the default parameters
    \param options (In/Out) : integer array, of size \c AZ_OPTIONS_SIZE, that will be
    populated with suitable values. A pointer to \c options will be stick into
    the parameters list. Note that this array is still required to apply the
    preconditioner! Do not delete options, nor let it go out of scope. The default value is
    0, meaning that \c SetDefaults() will allocate the array.
    \param params (In/Out) : double array, of size \c AZ_PARAMS_SIZE. See comments
    for \c options.
    \param OverWrite (In) : boolean.  If false, any pre-existing values in the
    parameter list will be preserved.  Default value is true, i.e., any
    pre-existing values may be overwritten.
   */
  int SetDefaults(std::string ProblemType, Teuchos::ParameterList & List,
		  int * options = 0, double * params = 0, const bool OverWrite=true);

  //! Sets default parameters for aggregation-based 2-level domain decomposition preconditioners.
  int SetDefaultsDD(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets default parameters for aggregation-based 2-level domain decomposition preconditioners, using LU on each subdomain
  int SetDefaultsDD_LU(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets default parameters for aggregation-based 3-level domain decomposition preconditioners.
  int SetDefaultsDD_3Levels(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets default parameters for aggregation-based 3-level domain decomposition preconditioners with LU.
  int SetDefaultsDD_3Levels_LU(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets default parameters for the eddy current equations equations.
  int SetDefaultsMaxwell(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets default parameters for classical smoothed aggregation.
  int SetDefaultsSA(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Sets defaults for energy minimization preconditioning for nonsymmetric problems.
  int SetDefaultsNSSA(Teuchos::ParameterList & List,
		               Teuchos::RCP<std::vector<int> > &options,
                       Teuchos::RCP<std::vector<double> > &params,
                       bool Overwrite=true);

  //! Reads in parameter list options from file.
  int ReadXML(const std::string &FileName, Teuchos::ParameterList &List,
                   const Epetra_Comm &Comm);

  //! Enumerated type indicating the type of AMG solver to be used.

  enum AMGType {

  ML_SA_FAMILY,         /*< Smoothed aggregation solver including EMIN,NSA */
  ML_MAXWELL,           /*< Old Maxwell solver */
  ML_COMPOSITE          /*< Composite AMG: block diagonal prolongator */
  };

/*!

   \brief MultiLevelPreconditioner: a class to define black-box multilevel preconditioners using aggregation methods.

   Class ML_Epetra::MultiLevelPreconditioner defined black-box algebraic
   multilevel preconditioners of matrices defined as Epetra_RowMatrix derived
   objects. The resulting preconditioner can be used in AztecOO, and in any
   other solver that accepts Epetra_Operator derived objects, and apply the
   action of the given Epetra_Operator using ApplyInverse().

   Please refer to the user's guide for a detailed introduction to
   this class, examples, and description of input parameters.

    This file requires ML to be configured with the following options:
    - \c --enable-epetra
    - \c --enable-teuchos

    The following option is suggested:
    - \c --enable-amesos
    - \c --enable-ifpack

    Some part of this class needs the following options:
    - \c --enable-aztecoo
    - \c --enable-anasazi

    It is important to note that ML is more restrictive than Epetra for
    the definition of maps. It is required that RowMatrixRowMap() is equal
    to OperatorRangeMap(). This is because ML needs to perform matrix-std::vector
    product, as well as getrow() functions, on the same data distribution.

    Also, for square matrices, OperatorDomainMap() must be as
    OperatorRangeMap().

    Several examples are provided in the \c examples subdirectories:
    - \ref ml_preconditioner_cpp is an introductory
      example;
    - \ref ml_2level_DD_cpp shows how to
      define a 2-level domain decomposition preconditioner using
      this class;
    - \ref ml_viz_cpp details how to visualize the aggregates;
    - \ref ml_maxwell_cpp reports how to
      use this class for Maxwell problems.

   \note
   Namespace ML_Epetra contains another Epetra_Operator derived class,
   ML_Epetra::MultiLevelOperator.
   - you should use MultiLevelOperator
     when your code already defines the required ML objects, with the optimal
     choice of parameters, and you just want to wrap the already defined ML
     preconditioners for AztecOO problems;
   - you should use MultiLevelPreconditioner
     when you have an Epetra_RowMatrix, and you don't want to code the
     conversion to ML_Operator, the creation of the hierarchy and the
     aggregates, and/or you want to experiment various combinations of the
     parameters, simply changing some parameters in a Teuchos::ParameterList.

   Defaults parameters can be specified using function SetDefaults().

    \author Marzio Sala, SNL 9214
*/
class MultiLevelPreconditioner : public virtual Epetra_Operator {

public:

  //@{ \name Constructors.

  //! Constructs a MultiLevelPreconditioner with default values.

  MultiLevelPreconditioner(const Epetra_RowMatrix & RowMatrix,
                           const bool ComputePrec = true);

  //! Constructs a MultiLevelPreconditioner. Retrieves parameters from \c List.

  MultiLevelPreconditioner(const Epetra_RowMatrix & RowMatrix,
			   const Teuchos::ParameterList & List,
			   const bool ComputePrec = true);

  //! Constructs a MultiLevelPreconditioner from an ML_Operator. Retrieves parameters from \c List.

  MultiLevelPreconditioner(ML_Operator* Operator,
			   const Teuchos::ParameterList& List,
			   const bool ComputePrec = true);

  //! Constructs a MultiLevelPreconditioner which is actually a composite AMG hierarchy using an array of ML_Operator's and an array of parameter lists.

  MultiLevelPreconditioner(ML_Operator *Operator,
                           const Teuchos::ParameterList& List,
                           Epetra_RowMatrix **DiagOperators,
			   Teuchos::ParameterList *DiagLists,
                           int NBlocks = 1,
			   const bool ComputePrec = true);

  //! \brief MultiLevelPreconditioner constructor for Maxwell's equations.
  /*! Takes the stiffness and mass terms of the matrix combined.

      \param EdgeMatrix - (In) Linear matrix to be solved.
      \param GradMatrix - (In) Node-to-edge connectivity matrix, a.k.a,
                               topological gradient
      \param NodeMatrix - (In) Auxiliary nodal finite element matrix
      \param List - (In) Teuchos parameter list containing solver options.
      \param ComputePrec - (In) Optional argument that specifies whether to
                                create preconditioner immediately.
                                Default is true.
      \param UseNodeMatrixForSmoother - (In) Use the nodal matrix for the nodal
                            portion of the Hipmair smoother (if used).
  */

  MultiLevelPreconditioner(const Epetra_RowMatrix& EdgeMatrix,
			   const Epetra_RowMatrix& GradMatrix,
			   const Epetra_RowMatrix& NodeMatrix,
			   const Teuchos::ParameterList& List,
			   const bool ComputePrec = true,
                           const bool UseNodeMatrixForSmoother = false);

  //! \brief MultiLevelPreconditioner constructor for Maxwell's equations.
  /*! Takes the stiffness and mass terms of the matrix separately.

      \param CurlCurlMatrix - (In) The curl-curl (stiffness) term of the
                                   matrix to be solved.
      \param MassMatrix - (In) The mass term of the matrix to be solved.
      \param GradMatrix - (In) Node-to-edge connectivity matrix, a.k.a,
                               topological gradient
      \param NodeMatrix - (In) Auxiliary nodal finite element matrix
      \param List - (In) Teuchos parameter list containing solver options.
      \param ComputePrec - (In) Optional argument that specifies whether to
                                create preconditioner immediately.
                                Default is true.
  */

  MultiLevelPreconditioner(const Epetra_RowMatrix & CurlCurlMatrix,
             const Epetra_RowMatrix & MassMatrix,
             const Epetra_RowMatrix & TMatrix,
             const Epetra_RowMatrix & NodeMatrix,
             const Teuchos::ParameterList & List,
             const bool ComputePrec = true);

#define NewStuff
#ifdef NewStuff
  //! Constructs a MultiLevelPreconditioner for multiphysics with variable dofs per node

  MultiLevelPreconditioner(Epetra_RowMatrix & RowMatrix,
			   const Teuchos::ParameterList & List,
                           const int & nNodes,
                           const int & maxDofPerNode,
                           bool * dofPresent,
                           Epetra_Vector & Lhs, 
                           Epetra_Vector & Rhs, 
                           const bool  rhsAndsolProvided,
			   const bool ComputePrec = true);
// ================================================ ====== ==== ==== == =
/*! Constructor for scalar PDE problems based on applying AMG to the distance
 *  Laplacian operator when constructing grid transfers. The main unique
 *  feature is that there may be some dofs that correspond to the same node
 *  location. These shared dofs fall into two categories. If these dofs are
 *  strongly connected to each other (as determined by tol), they are
 *  explicitly elminated from the Laplacian (merged into a supernode). Once
 *  a P is obtained, this P is then expanded to account for shared nodes
 *  by simply duplicating the supernodes row of P for each of the individual
 *  vertices that contribute to the supernode. If share dofs are weakly
 *  connected (or not connected at all), nothing special is done (other than
 *  the ususal ignoring of weak connections). One last trick is employed, 
 *  connections between supernodes and non-supernodes (i.e., regular nodes)
 *  are always assumed to be weak. Shared nodes are often used to capture
 *  interfaces or other features. By breaking these connections, the AMG
 *  can better maintain these features throughout the hierarchy. Finally, the
 *  constructor also allows for *  the removal of column nonzeros associated
 *  with Dirichlet points. To use this option the rhs and initial guess must be
 *  provided. Modification of the matrix, rhs, and initial guess must be 
 *  allowable to use this option.
 */

  MultiLevelPreconditioner(Epetra_RowMatrix & RowMatrix,
    const Teuchos::ParameterList & List,
    const double distTol, // two points are at the same location when
                          // || (x_1,y_1,z_1) -  (x_2,y_2,z_2)||_2 < distTol
    const double tol,     // ignore values when
                          //       A(i,j)^2 < A(i,i)*A(j,j)*tol^2
    Epetra_Vector & Lhs,
    Epetra_Vector & Rhs,
    const bool  rhsAndsolProvided,
    const bool ComputePrec = true);
#endif
#ifdef HAVE_ML_AZTECOO
  //! MultiLevelPreconditioner constructor for Maxwell's equations.
  /*! Takes the stiffness and mass terms of the matrix combined.  The edge
      matrix is of type Epetra_Msr, a light-weight wrapper for old-style Aztec
      MSR matrices.  This is intended as transition code for Aztec users.

      \param EdgeMatrix - (In) Linear matrix to be solved.
      \param GradMatrix - (In) Node-to-edge connectivity matrix, a.k.a,
                               topological gradient
      \param NodeMatrix - (In) Auxiliary nodal finite element matrix
      \param proc_config - (In) Aztec array specifying processor layout.
      \param List - (In) Teuchos parameter list containing solver options.
      \param ComputePrec - (In) Optional argument that specifies whether to
                                create preconditioner immediately.
                                Default is true.
  */

  MultiLevelPreconditioner(const Epetra_MsrMatrix & EdgeMatrix,
                         ML_Operator * GradMatrix,
                         AZ_MATRIX * NodeMatrix,
                         int       * proc_config,
                         const Teuchos::ParameterList & List,
                         const bool ComputePrec = true);
#endif

  //@}

  //@{ \name Destructor.

  //! Destroys the preconditioner.
  virtual ~MultiLevelPreconditioner() {
    if (IsComputePreconditionerOK_)
      DestroyPreconditioner();
  }

  //@}

  //@{ \name Query functions

  //! Prints label associated to this object.
  const char* Label() const{return(Label_);};

  //! Prints unused parameters in the input ParameterList on standard output.
  void PrintUnused() const
  {
    List_.unused(std::cout);
  }

  //! Prints unused parameters in the input ParameterList on the specified stream.
  void PrintUnused(std::ostream & os) const
  {
    List_.unused(os);
  }

  //! Prints unused parameters in the input ParameterList to std::cout on proc \c MyPID.
  /*! Mispelled parameters are simply ignored. Therefore, it is often the best
   * choice to print out the parameters that have not been used in the
   * construction phase.
   * - \param MyPID (In) : ID of process that should print the unused parameters.
   */
  void PrintUnused(const int MyPID) const;

  //! Gets a reference to the internally stored parameters' list.
  Teuchos::ParameterList& GetList()
  {
    return List_;
  }

  // Get a copy of the internally stored output list.
  Teuchos::ParameterList GetOutputList()
  {
    return OutputList_;
  }

  //! Prints on \c std::cout the values of the internally stored parameter list
  void PrintList();

  //! Copies \c List into the internally stored parameter list object.
  int SetParameterList(const Teuchos::ParameterList& List);

  //@}

  //@{ \name Mathematical functions.

  //! Apply the inverse of the preconditioner to an Epetra_MultiVector (NOT AVAILABLE)
  int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    return(-1);}

  //! Apply the preconditioner to an Epetra_MultiVector X, puts the result in Y
  int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

  //@}

  //@{ \name Attribute access functions


  //! Computes the multilevel hierarchy.
  /*! Computes the multilevel hierarchy. This function retrives the user's defines parameters (as
    specified in the input ParameterList), or takes default values otherwise, and creates the ML
    objects for aggregation and hierarchy. Allocated data can be freed used DestroyPreconditioner(),
    or by the destructor,

    In a Newton-type procedure, several linear systems have to be solved, Often, these systems
    are not too different. In this case, it might be convenient to keep the already
    computed preconditioner (with hierarchy, coarse solver, smoothers), and use it to
    precondition the next linear system. ML offers a way to determine whether the
    already available preconditioner is "good enough" for the next linear system.
    The user should proceed as follows:
    - define \c "reuse: enable" == \c true
    - solve the first linear system. ML tries to estimate the rate of convergence, and record it;
    - change the values of the linear system matrix (but NOT its structure)
    - compute the new preconditioner as \c ComputePreconditioner(true)
    It is supposed that the pointer to the Epetra_RowMatrix remains constant. Currently,
    it is not possible to modify this pointer (other than creating a new preconditioner)
  */

  int ComputePreconditioner(const bool CheckFiltering = false);

  /*! @brief Recompute the preconditioner (not implemented for Maxwell).

    @param[in] keepFineLevelSmoother : If true, the fine level smoother is not recomputed.  This is useful
    if the smoother is expensive to create, e.g., an incomplete factorization, and the fine level matrix
    has not changed.
  */
  int ReComputePreconditioner(bool keepFineLevelSmoother=false);

  //! Print the individual operators in the multigrid hierarchy.
  void Print(int level = -2);

  int ComputeAdaptivePreconditioner(int TentativeNullSpaceSize,
                                    double* TentativeNullSpace);

  //! Queries whether multilevel hierarchy has been computed or not.
  int IsPreconditionerComputed()  const
  {
    return(IsComputePreconditionerOK_);
  }

  // following functions are required to derive Epetra_RowMatrix objects.

  //! Sets ownership.
  int SetOwnership(bool ownership){ ownership_ = ownership; return(-1);};

  //! Sets use transpose (not implemented).
  int SetUseTranspose(bool useTranspose){return(-1);}

  //! Returns the infinity norm (not implemented).
  double NormInf() const {return(0.0);};

  //! Returns the current UseTranspose setting.
  bool UseTranspose() const {return(false);};

  //! Returns true if the \e this object can provide an approximate Inf-norm, false otherwise.
  bool HasNormInf() const{return(false);};

  //! Returns a pointer to the Epetra_Comm communicator associated with this operator.
  const Epetra_Comm& Comm() const{return(*Comm_);};

  //! Returns the Epetra_Map object associated with the domain of this operator.
  const Epetra_Map& OperatorDomainMap() const {return(*DomainMap_);};

  //! Returns the Epetra_Map object associated with the range of this operator.
  const Epetra_Map& OperatorRangeMap() const {return(*RangeMap_);};
  //@}

  //! Destroys all structures allocated in \c ComputePreconditioner() if the preconditioner has been computed.
  int DestroyPreconditioner();

  //! Returns a reference to the internally stored RowMatrix.
  const Epetra_RowMatrix& RowMatrix() const
  {
    return(*RowMatrix_);
  }

  //! Returns a reference to RowMatrix->Map().
  const Epetra_BlockMap& Map() const
  {
    return(RowMatrix_->Map());
  }

  //! Returns the global number of rows in the matrix.
  int NumGlobalRows() const
  {
    return(RowMatrix_->NumGlobalRows());
  }

  //! Returns the global number of columns in the matrix.
  int NumGlobalCols() const
  {
    return(RowMatrix_->NumGlobalCols());
  }

  //! Returns the local number of rows in the matrix.
  int NumMyRows() const
  {
    return(RowMatrix_->NumMyRows());
  }

  //! Returns the local number of columns in the matrix.
  int NumMyCols() const
  {
    return(RowMatrix_->NumMyCols());
  }

  //! Prints the computational stencil for the specified row and equation (for 2D Cartesian grids only)
  /*! For problems defined on 2D Cartesian grids (with node numbering increasing
   * along the x-axis), this function prints out the stencil in an intelligible
   * form.
   * \param nx (In) : number of nodes along the X-axis
   * \param ny (In) : number of nodes along the Y-axis
   * \param NodeID (In) : (local) ID of node that will be used to print the
   *   stencil. If set to -1, the code will automatically chose an internal node.
   *   Default: -1.
   * \param EquationID (In) : ID of the equation that will be used to print the
   *   stencil (default = 0)
   */
  int PrintStencil2D(const int nx, const int ny,
		     int NodeID = -1,
		     const int EquationID = 0);

  //! Cheap analysis of each level matrix.
  int AnalyzeHierarchy(const bool AnalyzeMatrices,
                       const int PreCycles, const int PostCycles,
                       const int MLCycles);

  //! Analyze the effect of each level's smoother on a random std::vector.
  int AnalyzeSmoothers(const int NumPreCycles = 1,
                       const int NumPostCycles = 1);

  //! Analyze the effect of the coarse solver on a random std::vector.
  int AnalyzeCoarse();

  //! Analyze the effect of the ML cycle on a random std::vector.
  int AnalyzeCycle(const int NumCycles = 1);

  //! Test several smoothers on fine-level matrix.
  int TestSmoothers(Teuchos::ParameterList& InputList,
		    const bool IsSymmetric = false);

  //! Test several smoothers on fine-level matrix using the current parameters.
  int TestSmoothers(const bool IsSymmetric = false) {
    return(TestSmoothers(List_,IsSymmetric));
  }

  //! Returns a pointer to the internally stored ml pointer
  const ML* GetML(const int WhichML= -1) const
  {
    if (WhichML < 0)
      return ml_;
    else if (WhichML == 0)
      return ml_nodes_;
    else
      return(0);
  }

  //! Returns a pointer to the internally stored agg pointer
  const ML_Aggregate* GetML_Aggregate() const
  {
    return agg_;
  }

  //! Generic interface to visualization methods.
  int Visualize(bool VizAggre, bool VizPreSmoother,
		bool VizPostSmoother, bool VizCycle,
		int NumApplPreSmoother, int NumApplPostSmoother,
		int NumCycleSmoother);

  //! Visualizes the shape of the aggregates.
  int VisualizeAggregates();

  //! Visualizes the effect of smoothers on a random std::vector.
  int VisualizeSmoothers(int NumPrecCycles = 1,
			 int NumPostCycles = 1);

  //! Visualizes the effect of the ML cycle on a random std::vector.
  int VisualizeCycle(int NumCycles = 1);

  /*! Creates label for this object (printed out by AztecOO).  This does not
      allocate/reallocate any memory.
  */
  int CreateLabel();

  void ReportTime();

  //! Return operator complexity and #nonzeros in fine grid matrix.
  void Complexities(double &complexity, double &fineNnz);

//@}

private:

  //! Copy constructor (NOT DEFINED)
  MultiLevelPreconditioner(const MultiLevelPreconditioner & rhs)
  {};

  //! operator = (NOT DEFINED)
  MultiLevelPreconditioner & operator = (const MultiLevelPreconditioner & rhs)
  {
    return *this;
  };

  //@{ \name Internal setting functions
  //! Initializes object with defauls values.
  int Initialize();

  /*! Sets smoothers.
    @param[in] skipFineLevelSmoother : If true, the fine level smoother is not set.  This is intended to be used in
    combination with ReComputePreconditioner.
  */
  int SetSmoothers(bool skipFineLevelSmoother=false);

  //! Sets coarse level solvers.
  int SetCoarse();

  //! Sets aggregation schemes.
  int SetAggregation();

  //! Sets preconditioner type (usually, V-cycle).
  int SetPreconditioner();

  //! Sets the null space for non-Maxwell problems.
  int SetNullSpace();

  //! Checks correctness of null space (discrete gradient) for Maxwell problems.
  //! The curl-curl and mass matrices must be supplied separately.
  void CheckNullSpace();

  //! Applies boundary conditions to gradient matrix.  (Maxwell's equations)
  void Apply_BCsToGradient( const Epetra_RowMatrix & EdgeMatrix,
                            const Epetra_RowMatrix & T);

  //! Transforms Epetra matrix column map (if necessary) to be compatible with
  /*! how ML handles column indices.  Any matrix that cannot be dynamically
      cast to an Epetra_CrsMatrix will not be changed.

      \param A - (In) Matrix that is to be transformed.
      \param transform - (In) EpetraExt widget that does the transformation.
      \param matrixName - (In) Optional label for the incoming matrix.
   */

#ifdef HAVE_ML_EPETRAEXT
  Epetra_RowMatrix* ModifyEpetraMatrixColMap( const Epetra_RowMatrix &A,
                                   EpetraExt::CrsMatrix_SolverMap &transform,
                                   const char* matrixName );
#endif
  //! Destroys Preconditioner if it not needed anymore. This includes some 'filtering' checks.

  int ConditionallyDestroyPreconditioner(const bool CheckPreconditioner);

  //! Set the finest level matrix in the MG hierarchy

  int SetFinestLevelMatrix();

  //! Set pointers indicating correspondence between array entries and MG levels

  int SetLevelIds(int Direction);

  //! Set eigenvalue scheme to be used by ML for spectral radius

  int SetEigenScheme();

  //! Dump various output matrices for debugging

  int MatrixDumper();

  //! Recompute complexities and print them.

  int ComputeAndPrintComplexities();

  //! Sets prolongator smoother parameters.
  int SetSmoothingDamping();

  //! Sets damping parameter for classical smoothed aggregation.
  int SetSmoothingDampingClassic();

#define OLD_AUX
#ifdef OLD_AUX
  int CreateAuxiliaryMatrixCrs(Epetra_FECrsMatrix * & FakeMatrix);

  int CreateAuxiliaryMatrixVbr(Epetra_VbrMatrix * & FakeMatrix);
#endif

  int SetupCoordinates();

  void PrintMem(char *fmt, int size, int, int);

  void PrintMemoryUsage();

  int SetFiltering();

  void RandomAndZero(double *, double *, int);

  //! Checks whether the previously computed preconditioner is still valuable for the newly available linear system.
  /*! Used only when \c "reuse: enable" is \c true, and
   * ComputePreconditioner(true) is called. */
  bool CheckPreconditionerKrylov();

  void VectorNorms(double*, int, double*,double*);

  //@}

  //@{ \name Internal data

  //! Pointer to ML_Struct
  ML* ml_;
  //! ML communicator, convenient to have separately from ml_,
  //  ml_nodes_, one or all of which may be null.
  ML_Comm* ml_comm_;

  //! indicates the type of AMG solver to be used: ML_SA_FAMILY, ML_MAXWELL, ML_COMPOSITE
  AMGType AMGSolver_;

  //! ML_Aggregate, contains aggregate information
  ML_Aggregate* agg_;
  //! Label for this object
  char* Label_;
  //! User-provided label for identifying preconditioner ctor/dtor, in the case
  //  of multiple instances of ML_Epetra::MultiLevelPreconditioner.
  std::string mlpLabel_;

  //! pointer to linear system matrix
  const Epetra_RowMatrix* RowMatrix_;

  //! AfineML_ points to the original ML operator passed in to the block
  //  matrix/composite version of the constructor.
  ML_Operator *AfineML_;

  //! Multigrid hierarchies applied to submatrices and used in a composite
  //  form to define the overall AMG hierarchy

  ML_Epetra::MultiLevelPreconditioner **SubMatMLPrec_;

  //! specifies whether a hierarchy already exists or not.
  bool IsComputePreconditionerOK_;

  //! Number of levels
  int NumLevels_;
  //! Domain Map
  const Epetra_Map* DomainMap_;
  //! Range Map
  const Epetra_Map* RangeMap_;
  //! Epetra communicator object
  const Epetra_Comm* Comm_;
  bool  ownership_;
  //! proc_config for Aztec smoothers
  int   ProcConfig_[AZ_PROC_SIZE];
  //! options for Aztec smoothers
  Teuchos::RCP<std::vector<int> >    SmootherOptions_;
  //! params for Aztec smoothers
  Teuchos::RCP<std::vector<double> > SmootherParams_;
  //! status for Aztec smoothers
  double SmootherStatus_[AZ_STATUS_SIZE];

  //! List containing all input parameters.
  Teuchos::ParameterList List_;
  //! List containing all output parameters
  Teuchos::ParameterList OutputList_;

  //! Maximum number of levels
  int MaxLevels_;

  //! Number of applications of the ML cycle
  int CycleApplications_;

  //! If \c true, zero starting solution is used in the application of the cycle.
  bool ZeroStartingSolution_;

  //! Integer array used to easily handle ML_INCREASING and ML_DECREASING
  /*! Integer array, of size MaxLevels_, that contain the ML level ID
    for the first logical level, and so on for all levels. The ML level ID
    of logical level L is LevelID_[L].
    In this interface, all levels move from 0 to MaxLevels-1.
    ML's level for interface's level i is LevelID_[i]
  */
  std::vector<int> LevelID_;

  //! If not NULL, contains the allocated null space std::vector
  double* NullSpaceToFree_;

  //! all std::cout's have this prefix (default'd in Initialize() )
  std::string PrintMsg_;
  //! all std::cerr's have this prefix (default'd in Initialize() )
  char ErrorMsg_[80];
  //! true if information has to be printed on this process
  bool verbose_;
  //! Number of PDE equations.
  int NumPDEEqns_;
  //! Number of iterations to use in profiling
  int profileIterations_;

  //@}

  //@{ \name Composite AMG variables

  //! Number of blocks making up composite operator

  int NBlocks_;

  //! Array of Diagonal Operators

  Epetra_RowMatrix **DiagOperators_;

  //! Array of Parameter lists for diagonal operators

  Teuchos::ParameterList *DiagLists_;

  //! special flag (currently for variable dof per node multiphysics AMG) skips smoother when building Laplacian hierarchy. 
  bool DontSetSmoothers_;
  //@}

  //@{ \name Maxwell variables

  //! Main matrix for Maxwell
  const Epetra_RowMatrix* EdgeMatrix_;
  //! stiffness and mass matrices
  const Epetra_RowMatrix* CurlCurlMatrix_;
  //! true if we summed curl-curl and mass
  bool CreatedEdgeMatrix_;
  const Epetra_RowMatrix* MassMatrix_;
  //! aux matrix for Maxwell
  const Epetra_RowMatrix* NodeMatrix_;
  //! T^T A T Matrix for use with Maxwell
  ML_Operator* TtATMatrixML_;
  bool UseNodeMatrixForSmoother_;
  bool CreatedNodeMatrix_;
  //! Auxiliary matrix used in intermediate step
  ML_Operator* ML_Kn_;
  bool CreatedML_Kn_;
  //! T matrix for Maxwell
  const Epetra_RowMatrix* TMatrix_;
#ifdef HAVE_ML_EPETRAEXT
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap RowMatrixColMapTrans_;
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap NodeMatrixColMapTrans_;
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap TMatrixColMapTrans_;
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap CurlCurlMatrixColMapTrans_;
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap MassMatrixColMapTrans_;
  //! Structure for compatibility between Epetra and ML column maps.
  EpetraExt::CrsMatrix_SolverMap TtATMatrixColMapTrans_;
#endif
  bool CreatedTMatrix_;
  ML_Operator* TMatrixML_;
  ML_Operator* TMatrixTransposeML_;
  ML_Operator** Tmat_array, ** Tmat_trans_array;
  ML_Operator** MassMatrix_array; // If curlcurl & mass are separate
  ML_Operator** CurlCurlMatrix_array;  // If curlcurl & mass are separate
  //! Auxiliary ML structure for Maxwell's equations.
  ML* ml_nodes_;

  void** nodal_args_,** edge_args_;

  //@}

  //@{ \name Variables for Timing
  //! Number of applications
  int NumApplications_;
  //! CPU time for all applications of the preconditioner
  double ApplicationTime_;
  bool FirstApplication_;
  //@ CPU time for first application
  double FirstApplicationTime_;
  //! Number of construction phases
  int NumConstructions_;
  //! CPU time for construction of the preconditioner.
  double ConstructionTime_;

  //@}

  // other stuff for old ML's compatibility
  Epetra_CrsMatrix* RowMatrixAllocated_;
  bool  AllocatedRowMatrix_; // used for composite constructor only

  bool AnalyzeMemory_;

  int memory_[ML_MEM_SIZE];

  // filtering stuff
  std::vector<double> flt_NullSpace_;
  ML* flt_ml_;
  ML_Aggregate* flt_agg_;

  // for reuse of preconditioning
  double RateOfConvergence_;

}; // class MultiLevelPreconditioner

} // namespace ML_Epetra

#ifdef NewStuff
/******************************************************************************
  Vector class that is used for for multiphysics multigrid on matrices with
  varaible dofs per node. Didn't want anything specific to epetra, ML, and 
  MueLu; didn't like new/malloc features of std::vector, and I have a hard
  time figuring out teuchos stuff ... so I implemented my own. The syntax
  is pretty similar to std::vector, but it is more efficient from a memory
  allocation perspective (though less flexible). 
******************************************************************************/

enum free_s {useFree, useDelete, neither};
template <class T> class MLVec {
  public:

    // different constructors & destructors

    MLVec();              
    MLVec(int mysize);
    MLVec(T *startptr, T *endptr, bool okToChangePtr = false, free_s freeType = neither);
    ~MLVec();

    int     size() const;        // return vector size
    T       *getptr();           // return pointer to start of raw vector data
    void    wrap(T *startptr, T *endptr,// same as above construct but allows
            bool okToChangePtr = false, // one to first create an empty MLVec
            free_s freeType = neither); // and then wrap it with data.
                                 // is not invoked when MLVec goes away
    void    relinquishData();    // disassociates data from MLVec, so free 
                                 // is not invoked when MLVec goes away
    void    resize(int newSize); // change vector size. If newSize is smaller
                                 // than size(), keep first newSize values
    bool    empty() const;       // indicates whether vector is empty

    inline T& operator [] (int index) { return data_[index]; }
    inline const T& operator [] (int index) const { return data_[index]; }

  private:
    int    size_;
    T      *data_;
    bool   okToChangePtr_;   // These last few indicate the status of vector.
    free_s freeType_;        // This is used to see if memory can be freed, 
    bool   getPtrInvoked_;   // reallocated, or pointers can be moved. 
};

// Create an empty vector 

template <class T> MLVec<T>::MLVec() {
  data_          = NULL;      size_          = 0;
  okToChangePtr_ = true;      getPtrInvoked_ = false;
  freeType_      = neither;
}

// Create an vector without any data but allocate space
template <class T> MLVec<T>::MLVec(int mysize) {

  TEUCHOS_TEST_FOR_EXCEPTION(mysize < 0,std::logic_error,
  "MLVec error, cannot create a negative size (= " << mysize << ") vector\n");

  data_          = NULL;        size_          = mysize;
  okToChangePtr_ = true;        getPtrInvoked_ = false;
  freeType_      = neither;

  if (size_ > 0) {
     data_ = (T *) ML_allocate(sizeof(T)*size_); 

     TEUCHOS_TEST_FOR_EXCEPTION(data_ == NULL,std::logic_error,
     "MLVec error, not enough space for vector of size " << size_ << "\n");

     freeType_ = useFree;
  }
}

// Create a vector that wraps existing data that starts and ends at 'start'
// and 'end' respectively. In most cases, this means that the pointers
// associated with the raw data cannot be changed, and that the raw data
// should remain when this vector goes away 
template <class T> MLVec<T>::MLVec(T *start, T *end, bool okToChangePtr,
                                   free_s freeType) { 

  TEUCHOS_TEST_FOR_EXCEPTION(end < start,std::logic_error,
  "MLVec error, cannot create a negative size (= " << end-start<< ") vector\n");

  freeType_      = freeType;   okToChangePtr_ = okToChangePtr;
  data_          = start;      size_          = end - start;
  if (size_ == 0) data_ = NULL;
  getPtrInvoked_ = false;
}
template <class T> void MLVec<T>::wrap(T *start, T *end, bool okToChangePtr,
                                   free_s freeType) { 

  TEUCHOS_TEST_FOR_EXCEPTION((data_ != NULL)||(size_  != 0)||(okToChangePtr_ != true)||
                             (getPtrInvoked_ != false)||(freeType_!= neither),
                             std::logic_error,
  "MLVec wrap error, wrap() must be invoked on vectors created with default constructor\n");

  TEUCHOS_TEST_FOR_EXCEPTION(end < start,std::logic_error,
  "MLVec error, cannot create a negative size (= " << end-start<< ") vector\n");

  freeType_      = freeType;   okToChangePtr_ = okToChangePtr;
  data_          = start;      size_          = end - start;
  if (size_ == 0) data_ = NULL;
  getPtrInvoked_ = false;
}

template <class T> MLVec<T>::~MLVec()              { 
   if (data_ != NULL) {
      TEUCHOS_TEST_FOR_EXCEPTION(size_ == 0,std::logic_error,
      "MLVec error, data pointer should be null for 0 length vectors\n");
      if (okToChangePtr_ && (freeType_ == useFree))   ML_free(data_); 
      if (okToChangePtr_ && (freeType_ == useDelete)) delete data_;
   }
}

// data_ is no longer associated with this vector. The raw data continues to
// exist.  Can be used in conjunction with getptr() to first get the raw
// data out of MLVec and then effectively empty MLVec.
template <class T> void MLVec<T>::relinquishData() { 

  TEUCHOS_TEST_FOR_EXCEPTION((data_ != NULL) && (getPtrInvoked_ == false),
  std::logic_error,"relinquishData without invoking getPtr() should lead"
  " to memory leak\n");

  data_          = NULL;     size_          = 0; 
  okToChangePtr_ = true;     getPtrInvoked_ = false;
  freeType_      = neither;
}
template <class T> int MLVec<T>::size() const  {   return size_; }
template <class T> T* MLVec<T>::getptr() {getPtrInvoked_ =true; return data_;}
template <class T> bool MLVec<T>::empty() const {if (size_ > 0) return false; 
                                                 else return true;}

// change the size of the vector. There are lots of restrictions as to when 
// this is allowed. There are also lots of cases depending on whether the 
// original vector was zero length (i.e. just allocate memory), the new vector
// is zero length (i.e. just free memory), the new vector is longer or shorter
// than the original vector.
template <class T> void MLVec<T>::resize(int newsize) {

  if (newsize == size_) return;

  TEUCHOS_TEST_FOR_EXCEPTION( newsize < 0, 
  std::logic_error,"MLVec error, cannot resize to a negative length (= " << 
  newsize << ") vector \n");

  TEUCHOS_TEST_FOR_EXCEPTION((newsize > size_) && (size_ > 0) && 
                             (freeType_ == useDelete), std::logic_error,
  "MLVec error, cannot resize vector allocated with new\n");

  if ( (newsize <= size_) && (size_ > 0) && (freeType_ == useDelete) &&
       (okToChangePtr_) ) { size_ = newsize; return; }

  if (okToChangePtr_) {
     if (size_ == 0) {
       if ( newsize > 0) {
         data_ = (T *) ML_allocate(sizeof(T)*newsize);

         TEUCHOS_TEST_FOR_EXCEPTION( data_ == NULL, std::logic_error,
         "MLVec error, not enough resize space for " << newsize <<
         " length vector\n");

         freeType_ = useFree;
       }
       else data_ = NULL;
     }
     else {
       if ( newsize > 0) {
          T* newStuff = (T *) ML_allocate(sizeof(T)*newsize);
          if (size_ < newsize) 
             for (int i = 0; i < size_; i++) newStuff[i] = data_[i];
          else 
             for (int i = 0; i < newsize; i++) newStuff[i] = data_[i];
          ML_free(data_);
          data_ = newStuff;
         TEUCHOS_TEST_FOR_EXCEPTION( data_ == NULL, std::logic_error,
         "MLVec error, not enough resize realloc space for " << newsize <<
         " length vector\n");

       }
       else {
          if (freeType_ == useFree) { ML_free(data_); data_ = NULL; }
          else if (freeType_ == useDelete) { delete data_; data_ = NULL; }

         TEUCHOS_TEST_FOR_EXCEPTION( freeType_ == neither, std::logic_error,
         "MLVec error, do not know how to free this vector\n"); 
        
       }
     }
     size_ = newsize;
  }
  TEUCHOS_TEST_FOR_EXCEPTION(okToChangePtr_ == false, std::logic_error,
  "MLVec error, not allowed to resize this type of wrapped vector\n"); 
}


/******************************************************************************
  Data structures & functions to hide communication details associated with 
  epetra matrices, ML matrices, or MueLu matrices (i.e.  xpetra operators). 
******************************************************************************/

enum epetraOrMLOrMueLu {epetraType , mlType, mueluType}; 

struct wrappedCommStruct {
    epetraOrMLOrMueLu whichOne;
    void                 *data;  // matrix on which communication is based
    int                 nProcs;  // total number of processors participating
    int                  myPid;  
    int          maxDofPerNode;
    int                vecSize;  
};

template <class T> int nodalComm(MLVec<T>& vector, MLVec<int>& myLocalNodeIds,
                                 struct wrappedCommStruct& framework)
{
   /***************************************************************************
    Performs communication to update an amalgamated/nodal vector using either
    an ML style or an epetra style unamalgamated/dof matrix to define the
    communication/import pattern. To do this, it utilzes myLocalNodeIds[i]
    which indicates that the ith dof lies within the myLocalNodeIds[i]th
    node. The ghost portion of myLocalNodeIds[] was computed by 
    assignGhostLocalNodeIds().

    Note: As several dofs lie within each node, several of myLocalNodeIds[]'s
    entries will be duplicates. 
   ***************************************************************************/

    MLVec<T> temp(myLocalNodeIds.size());

    // copy from nodal vector to dof vector
    for (int i = 0; i < myLocalNodeIds.size(); i++)
       temp[i] = vector[ myLocalNodeIds[i]];

    dofComm(temp, framework);

    // copy from dof vector to nodal vector
    for (int i = 0; i < myLocalNodeIds.size(); i++)
       vector[ myLocalNodeIds[i]] = temp[i];

    return 0;
}

template <class T> int dofComm(MLVec<T>& vector,
                               struct wrappedCommStruct& framework)
{
   /***************************************************************************
    Performs communication to update vector using either an ML style matrix
    or an epetra style matrix to define the communication/import pattern.

    Note: Directly utilizes the communication of the underlying matrix
    as opposed to nodalComm() and dofCommUsingMlNodalMatrix(). nodalComm() uses
    an unamalgamated matrix to define the communication associated with an
    amalgamated matrix. dofCommUsingMlNodalMatrix() uses an amalgamated/nodal
    matrix maxDofPerNode times to interpolate each dof of an unamalgamated
    matrix. dofCommUsingMlNodalMatrix() assumes all nodes have maxDofPerNode
    degrees-of-freedom (so it is only appropriate on coarser levels).
   ***************************************************************************/

   MLVec<double> aCopy(vector.size());

   // need to make a copy of the data as we will only use framework comm()
   // functions that work with doubles

   for (int i = 0; i < vector.size(); i++) aCopy[i] = (double) vector[i];

   double *vectorData = aCopy.getptr();

   if (framework.whichOne == mlType) {

      ML_Operator *Amat = (ML_Operator *) framework.data;

      // grab the data pointer

      if (Amat == NULL) return 0;
      if (Amat->getrow == NULL) return 0;
      if (Amat->getrow->pre_comm == NULL) return 0;

      ML_exchange_bdry(vectorData,Amat->getrow->pre_comm,
                 Amat->invec_leng,Amat->comm,ML_OVERWRITE,NULL);
   }
   else {
      Epetra_CrsMatrix *Amat = (Epetra_CrsMatrix *) framework.data;
      Epetra_RowMatrix *ArMat = (Epetra_RowMatrix *) Amat;

      ML_Epetra_comm_wrapper(vectorData, (void *) ArMat);
   }

   for (int i = 0; i < vector.size(); i++) vector[i] = (T) aCopy[i];

   return(0);
}

extern "C" int dofCommUsingMlNodalMatrix(double *data, void *widget);

extern int MLsortCols(const MLVec<int>& ARowPtr, MLVec<int>& ACols, 
                    MLVec<double>& AVals);

extern int MLnMyGhost(struct wrappedCommStruct& framework);

extern int MLfillNodalMaps(MLVec<int> &amalgRowMap, MLVec<int> &amalgColMap,
        MLVec<int> &myLocalNodeIds, int nLocalDofs, 
        struct wrappedCommStruct &framework,
        int nLocalNodes, int nLocalPlusGhostNodes);

extern int MLcolGlobalIds(struct wrappedCommStruct& framework, MLVec<int>& myGids);

extern int MLassignGhostLocalNodeIds(MLVec<int> &myLocalNodeIds, int nLocalDofs,
                   int nLocalPlusGhostDofs, struct wrappedCommStruct &framework,
                   int &nLocalNodes, int &nLocalPlusGhostNodes);

extern int MLextractDiag(const MLVec<int>& rowPtr, const MLVec<int>& cols,
                       const MLVec<double>& , MLVec<double>& diagonal,
                       struct wrappedCommStruct &framework);

extern int MLfindDirichlets(const MLVec<int>& rowPtr, const MLVec<int>& cols,
                          const MLVec<double>& vals,
                          const MLVec<double>& diagonal,
                          double tol, MLVec<bool>& dirOrNot,
                          struct wrappedCommStruct &framework);

extern int MLrmDirichletCols(MLVec<int>& rowPtr, MLVec<int>& cols,
                           MLVec<double>& vals, const MLVec<double>& diagonal,
                           bool squeeze, MLVec<double>& solution,
                           MLVec<double>& rhs, const MLVec<bool>& dirOrNot, 
                           struct wrappedCommStruct &framework);

extern int MLsqueezeOutNnzs(MLVec<int>& rowPtr, MLVec<int>& cols,
                          MLVec<double>& vals, const MLVec<bool>& keep);

extern int MLbuildMap(const MLVec<bool>& dofPresent, MLVec<int>& map, int NDof);

extern int MLvariableDofAmalg(int nCols, const MLVec<int>& rowPtr,
                          const MLVec<int>& cols, const MLVec<double>& vals,
                          int nNodes, int maxDofPerNode, const MLVec<int>& map,
                          const MLVec<double>& diag, double tol,
                          MLVec<int>& amalgRowPtr, MLVec<int>& amalgCols,
                          struct wrappedCommStruct &framework,
                          MLVec<int>& myLocalNodeIds);


extern int MLrmDifferentDofsCrossings(const MLVec<bool>& dofPresent,
                                    int maxDofPerNode, MLVec<int>& rowPtr,
                                    MLVec<int>& cols, int nCols,
                                    struct wrappedCommStruct& framework, 
                                    MLVec<int> &myLocalNodeIds);

extern int MLbuildLaplacian(const MLVec<int>& rowPtr, const MLVec<int>& cols,
                          MLVec<double>& vals, const MLVec<double>& x,
                          const MLVec<double>& y, const MLVec<double>& z);

extern int MLunamalgP(const MLVec<int>& amalgRowPtr,
                    const MLVec<int>& amalgCols,
                    const MLVec<double>& amalgVals, int maxDofPerNode,
                    const MLVec<char>& status, bool fineIsPadded,
                  MLVec<int>& rowPtr, MLVec<int>& cols, MLVec<double>& vals);

extern int MLfindEmptyRows(const MLVec<int>& rowPtr, const MLVec<int>& cols,
                         const MLVec<double>& vals, MLVec<bool>& rowEmpty);

extern int MLreplaceEmptyByDirichlet(MLVec<int>& rowPtr, MLVec<int>& cols,
                        MLVec<double>& vals, const MLVec<bool>& colEmpty);

extern int MLfineStatus(const MLVec<bool>& dofPresent,
                           const MLVec<int>& map, const MLVec<bool>& dirOrNot,
                           MLVec<char>& status);

extern int MLcoarseStatus(const MLVec<bool>& rowEmpty,
                        const MLVec<bool>& dirOrNot, MLVec<char>& status);



extern int MLShove(ML_Operator *Mat, MLVec<int>& rowPtr, MLVec<int>& cols, MLVec<double>& vals, int invec_leng, int (*commfunc  )(double *vec, void *data), struct wrappedCommStruct& framework, int nGhost);


extern int ZeroDist(MLVec<double>& xxx, MLVec<double>& yyy, MLVec<double>& zzz,
             MLVec<int>& rowPtr, MLVec<int>& cols, MLVec<double>& vals,
             MLVec<double>& diagonal, double tol, MLVec<int>& rowZeros, MLVec<int>& colZeros,
             double disttol);

extern int MergeShared(MLVec<int>& cols, MLVec<int>& rowZeros, MLVec<int>& colZeros, 
                       MLVec<int>& groupHead, MLVec<int>& groupNext);

extern int BuildNonSharedMap(MLVec<int>& newMap, MLVec<int>& groupHead, MLVec<int>& groupNext);

extern int buildCompressedA(MLVec<int>& inputRowPtr, MLVec<int>& inputCols,
                     MLVec<double>& inputVals, MLVec<double>& diagonal,
                     MLVec<int>& groupHead, MLVec<int>& groupNext, 
                     MLVec<int>& outputRowptr, MLVec<int>& outputCols,
                     MLVec<int>& map, int newN, double tol);


#endif /* NewStuff */
#endif /* defined HAVE_ML_EPETRA and HAVE_ML_TEUCHOS */

#endif /* define ML_MULTILEVELPRECONDITIONER_H */