{ "cells": [ { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true, "slideshow": { "slide_type": "slide" } }, "source": [ "# MueLu - Level smoothers" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true, "slideshow": { "slide_type": "slide" } }, "source": [ "## General remarks" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true, "slideshow": { "slide_type": "slide" } }, "source": [ "- *MueLu* uses algorithms from other *Trilinos* packages for level smoothers. These include *Ifpack/Ifpack2* and *Teko*. Similar for direct solvers: *MueLu* uses the direct solver provided by *Amesos/Amesos2*.\n", "- For creating concrete smoothers on the different levels we use the *prototype* pattern. I.e. we define one (or more) smoother prototypes which are given to the *SmootherFactory* class responsible for the concrete level." ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Level smoothers in MueLu" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### General preparations" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Declare location of include headers " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ ".I /opt/install/do-conf-ep-serial/include/" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ " Load Trilinos libraries" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ ".L libepetra" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ ".L libxpetra" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ ".L libmuelu" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ ".L libteuchoscore" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Include some standard header files to create a communicator and an ```Xpetra::Map```" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"Kokkos_DefaultNode.hpp\"\n", "#include \"Epetra_SerialComm.h\"\n", "#include \"Epetra_Vector.h\"\n", "#include \"Epetra_CrsMatrix.h\"\n", "#include \"Teuchos_RCP.hpp\"\n", "#include \"Teuchos_DefaultSerialComm.hpp\"\n", "#include \"Xpetra_MapFactory.hpp\"\n", "#include \"Xpetra_Map.hpp\"\n", "#include \"Xpetra_EpetraCrsMatrix.hpp\"" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Again, we use standard template parameters:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "typedef double Scalar;\n", "typedef int LocalOrdinal;\n", "typedef int GlobalOrdinal;\n", "typedef Kokkos::Compat::KokkosSerialWrapperNode EpetraNode;\n", "typedef EpetraNode Node;\n", "typedef Scalar SC;\n", "typedef LocalOrdinal LO;\n", "typedef GlobalOrdinal GO;\n", "typedef Node NO;" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"Xpetra_CrsMatrixWrap.hpp\"\n", "#include \"Xpetra_CrsMatrix.hpp\"\n", "#include \"Xpetra_MultiVector.hpp\"\n", "#include \"Xpetra_BlockedMultiVector.hpp\"\n", "#include \"Xpetra_BlockedVector.hpp\"\n", "#include \"Xpetra_MultiVectorFactory.hpp\"" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "### Create a matrix object" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Let's create a matrix and wrap it as an ```Xpetra::CrsMatrixWrap``` object" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(int) 30\n" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "int NumMyElements = 30;" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(int) 0\n" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Epetra_SerialComm Comm;\n", "Epetra_Map Map(-1, NumMyElements, 0, Comm);\n", "int NumGlobalElements = Map.NumGlobalElements();\n", "\n", "Epetra_CrsMatrix A(Copy, Map, 3);\n", "\n", "double negOne = -1.0;\n", "double posTwo = 2.0;\n", "double posOne = +1.0;\n", "double zero = 0.0;\n", "for (int i=0; i > &) @0x7f4a182be318\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Teuchos::RCP rcpA = Teuchos::rcpFromRef(A);\n", "Xpetra::EpetraCrsMatrixT xA(rcpA);\n", "Teuchos::RCP > xCrsMat = Teuchos::rcpFromRef(xA);\n", "Xpetra::CrsMatrixWrap xCrsWrap = Xpetra::CrsMatrixWrap(xCrsMat);\n", "Teuchos::RCP > xMat = Teuchos::rcpFromRef(xCrsWrap);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "The ```xCrsWrap``` object contains the matrix object that can be used as input for *MueLu*." ] }, { "cell_type": "markdown", "metadata": { "collapsed": false, "deletable": true, "editable": true }, "source": [ "### Create a smoother" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "Smoother details are stored in a ```Teuchos::ParameterList```. For demonstration purposes we create a Symmetric Gauss-Seidel smoother ($50$ sweeps, damping factor $0.45$):" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"Teuchos_ParameterList.hpp\"" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::ParameterList &) {}\n" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Teuchos::ParameterList smooParams;\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::ParameterList &) { \"relaxation: type\" => @0xc1fda48, \"relaxation: sweeps\" => @0xc1fda90, \"relaxation: damping factor\" => @0xc1fdad8 }\n" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smooParams.set(\"type\", Teuchos::ParameterEntry(\"RELAXATION\"));\n", "Teuchos::ParameterList& subSmooParams = smooParams.sublist(\"ParamterList\");\n", "subSmooParams.set(\"relaxation: type\", \"Symmetric Gauss-Seidel\");\n", "subSmooParams.set(\"relaxation: sweeps\",50);\n", "subSmooParams.set(\"relaxation: damping factor\",0.45);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "You can easily print the content of a parameter list using" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "type = RELAXATION [unused] [unused]\n", "ParamterList -> \n", " relaxation: type = Symmetric Gauss-Seidel [unused]\n", " relaxation: sweeps = 50 [unused]\n", " relaxation: damping factor = 0.45 [unused]\n", "\n" ] }, { "data": { "text/plain": [ "(std::basic_ostream >::__ostream_type &) @0x7f4a10003700\n" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "std::cout << smooParams << std::endl;" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "The ```TrilinosSmoother``` class is a wrapper class for most relaxation-based and ILU type smoothers (provided by *Ifpack* or *Ifpack2*):" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_TrilinosSmoother.hpp\"" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "To create a new ```TrilinosSmoother``` you have to provide a type (e.g., ```RELAXATION```) and the corresponding *Ifpack/Ifpack2* parameter list:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::TrilinosSmoother &) @0x7f4a182be3e0\n" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::TrilinosSmoother trSmoother(\"RELAXATION\", subSmooParams);\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "The smoother use variable ```A``` stored on the ```Level``` class as input. Therefore, let's create a ```Level``` container and store the matrix:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_Level.hpp\"" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::Level &) @0x7f4a182be5c0\n" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::Level l;" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "
Note:It is important, that ```A``` is a ```RCP>``` object for the smoother to be able to properly cast the matrix object (and determine whether we are using the *Epetra* or *Tpetra* stack).\n", "
" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l.Set(\"A\",xMat);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Furthermore, we have to store the type of the underlying linear algebra package used underneath in the ```Level``` class. This is necessary such that the ```TrilinosSmoother``` class can distinguish between *Epetra* and *Tpetra* and create either *Ifpack* or *Ifpack2* preconditioners.\n", "
Note:We have to explicitly set this flag here. When using the standard MueLu setup routines, this is automatically done and the user does not have to think about that.\n", "
" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l.setlib(Xpetra::UseEpetra);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then, we can request the ```Smoother```. Without requesting the smoother the smoother would not be built." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l.Request(\"Smoother\", &trSmoother);" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "LevelID = -1\n", " data name gen. factory addr. req keep type data req'd by \n", " -------------------- ------------------ ----- ----- --------------- -------------- --------------------\n", " A NoFactory 1 User Matrix available 0x11f22808 \n", " Smoother 0x7f4a182be3e8 1 No unknown not available 0x123b69e0 \n" ] }, { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l.print(std::cout, MueLu::Extreme);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then we can call the ```Setup``` routine of the smoother object. Depending on the smoother type, the ILU triangle block are built or the diagonal is inverted. We only have to provide the current level container with ```A``` stored." ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Setup Smoother (MueLu::IfpackSmoother{type = point relaxation stand-alone})\n", " IFPACK (Local SGS, sweeps=50, damping=0.45)\n" ] }, { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trSmoother.Setup(l);" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "For using the smoother we need a solution guess and a RHS vector." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::RCP > &) @0x7f4a182be6c0\n" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Teuchos::RCP> x = Xpetra::MultiVectorFactory::Build(xMat->getRowMap(),1);" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::RCP > &) @0x7f4a182be6d8\n" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Teuchos::RCP> b = Xpetra::MultiVectorFactory::Build(xMat->getRowMap(),1);" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b->putScalar(1.0);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then we can call the ```Apply``` routine. In our example it performs $50$ SGS sweeps (damping $0.45$):" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trSmoother.Apply(*x,*b,true);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then you print the solution vector and the matrix..." ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " MyPID GID Value \n", " 0 0 1\n", " 0 1 6.49166\n", " 0 2 11.0806\n", " 0 3 14.8804\n", " 0 4 17.9917\n", " 0 5 20.51\n", " 0 6 22.5236\n", " 0 7 24.1134\n", " 0 8 25.3514\n", " 0 9 26.3004\n", " 0 10 27.0142\n", " 0 11 27.5373\n", " 0 12 27.9054\n", " 0 13 28.1456\n", " 0 14 28.2765\n", " 0 15 28.3089\n", " 0 16 28.2461\n", " 0 17 28.0834\n", " 0 18 27.8084\n", " 0 19 27.4012\n", " 0 20 26.8334\n", " 0 21 26.0683\n", " 0 22 25.0607\n", " 0 23 23.7564\n", " 0 24 22.0922\n", " 0 25 19.9963\n", " 0 26 17.3887\n", " 0 27 14.1817\n", " 0 28 10.2813\n", " 0 29 5.58836\n", "\n" ] }, { "data": { "text/plain": [ "(std::basic_ostream >::__ostream_type &) @0x7f4a10003700\n" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "std::cout << *x << std::endl;" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Epetra::CrsMatrix\n", "Number of Global Rows = 30\n", "Number of Global Cols = 30\n", "Number of Global Diagonals = 30\n", "Number of Global Nonzeros = 88\n", "Global Maximum Num Entries = 3\n", "\n", "\n", "\n", "\n", "\n", "Number of My Rows = 30\n", "Number of My Cols = 30\n", "Number of My Diagonals = 30\n", "Number of My Nonzeros = 88\n", "My Maximum Num Entries = 3\n", "\n", " Processor Row Index Col Index Value \n", " 0 0 0 1 \n", " 0 0 1 0 \n", " 0 1 0 -1 \n", " 0 1 1 2 \n", " 0 1 2 -1 \n", " 0 2 1 -1 \n", " 0 2 2 2 \n", " 0 2 3 -1 \n", " 0 3 2 -1 \n", " 0 3 3 2 \n", " 0 3 4 -1 \n", " 0 4 3 -1 \n", " 0 4 4 2 \n", " 0 4 5 -1 \n", " 0 5 4 -1 \n", " 0 5 5 2 \n", " 0 5 6 -1 \n", " 0 6 5 -1 \n", " 0 6 6 2 \n", " 0 6 7 -1 \n", " 0 7 6 -1 \n", " 0 7 7 2 \n", " 0 7 8 -1 \n", " 0 8 7 -1 \n", " 0 8 8 2 \n", " 0 8 9 -1 \n", " 0 9 8 -1 \n", " 0 9 9 2 \n", " 0 9 10 -1 \n", " 0 10 9 -1 \n", " 0 10 10 2 \n", " 0 10 11 -1 \n", " 0 11 10 -1 \n", " 0 11 11 2 \n", " 0 11 12 -1 \n", " 0 12 11 -1 \n", " 0 12 12 2 \n", " 0 12 13 -1 \n", " 0 13 12 -1 \n", " 0 13 13 2 \n", " 0 13 14 -1 \n", " 0 14 13 -1 \n", " 0 14 14 2 \n", " 0 14 15 -1 \n", " 0 15 14 -1 \n", " 0 15 15 2 \n", " 0 15 16 -1 \n", " 0 16 15 -1 \n", " 0 16 16 2 \n", " 0 16 17 -1 \n", " 0 17 16 -1 \n", " 0 17 17 2 \n", " 0 17 18 -1 \n", " 0 18 17 -1 \n", " 0 18 18 2 \n", " 0 18 19 -1 \n", " 0 19 18 -1 \n", " 0 19 19 2 \n", " 0 19 20 -1 \n", " 0 20 19 -1 \n", " 0 20 20 2 \n", " 0 20 21 -1 \n", " 0 21 20 -1 \n", " 0 21 21 2 \n", " 0 21 22 -1 \n", " 0 22 21 -1 \n", " 0 22 22 2 \n", " 0 22 23 -1 \n", " 0 23 22 -1 \n", " 0 23 23 2 \n", " 0 23 24 -1 \n", " 0 24 23 -1 \n", " 0 24 24 2 \n", " 0 24 25 -1 \n", " 0 25 24 -1 \n", " 0 25 25 2 \n", " 0 25 26 -1 \n", " 0 26 25 -1 \n", " 0 26 26 2 \n", " 0 26 27 -1 \n", " 0 27 26 -1 \n", " 0 27 27 2 \n", " 0 27 28 -1 \n", " 0 28 27 -1 \n", " 0 28 28 2 \n", " 0 28 29 -1 \n", " 0 29 28 -1 \n", " 0 29 29 2 \n", "\n" ] }, { "data": { "text/plain": [ "(std::basic_ostream >::__ostream_type &) @0x7f4a10003700\n" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "std::cout << A << std::endl;" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "... or you calculate the Residual or Residual norm:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"Teuchos_Array.hpp\"" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_Utilities.hpp\"" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::Array::Magnitude>) { 3.606524 }\n" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::Utilities::ResidualNorm(*xMat,*x,*b);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Create a direct solver object" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "In a similar way, we can define a direct solver. As type, we choose a direct solver from *Amesos/Amesos2* such as ```KLU``` and a parameter list. In our example an empty parameter list is sufficient." ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_DirectSolver.hpp\"" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::DirectSolver &) @0x7f4a182be790\n" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Teuchos::ParameterList dsPl;\n", "MueLu::DirectSolver ds(\"Klu\", dsPl);" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x->putScalar(0.0);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Before we call the ```Setup``` routine we have to request the ```DirectSolver``` on the level class." ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l.Request(\"DirectSolver\", &ds);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then, it is time for the ```Setup``` call which usually performs the factorization." ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Setup Smoother (MueLu::AmesosSmoother{type = Klu})\n" ] }, { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds.Setup(l);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Now, we can call the ```Apply``` routine and solve the problem (exactly)." ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds.Apply(*x,*b,true);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Again, we can print out the solution and calculate the residual norm:" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " MyPID GID Value \n", " 0 0 1\n", " 0 1 15.4667\n", " 0 2 28.9333\n", " 0 3 41.4\n", " 0 4 52.8667\n", " 0 5 63.3333\n", " 0 6 72.8\n", " 0 7 81.2667\n", " 0 8 88.7333\n", " 0 9 95.2\n", " 0 10 100.667\n", " 0 11 105.133\n", " 0 12 108.6\n", " 0 13 111.067\n", " 0 14 112.533\n", " 0 15 113\n", " 0 16 112.467\n", " 0 17 110.933\n", " 0 18 108.4\n", " 0 19 104.867\n", " 0 20 100.333\n", " 0 21 94.8\n", " 0 22 88.2667\n", " 0 23 80.7333\n", " 0 24 72.2\n", " 0 25 62.6667\n", " 0 26 52.1333\n", " 0 27 40.6\n", " 0 28 28.0667\n", " 0 29 14.5333\n", "\n" ] }, { "data": { "text/plain": [ "(std::basic_ostream >::__ostream_type &) @0x7f4a10003700\n" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "std::cout << *x << std::endl;" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(Teuchos::Array::Magnitude>) { 0.000000 }\n" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::Utilities::ResidualNorm(*xMat,*x,*b);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "## Level smoothers in a MueLu hierarchy" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "In order to make the smoothers available to the multigrid hierarchy, we have to create\n", "1. the smoothers\n", "2. Create a ```SmootherFactory``` class accepting the smoother objects for pre- and/or postsmoothing\n", "3. Register the ```SmootherFactory``` object in a ```FactoryManager```\n", "\n", "More details below and in the ```Hierarchy``` lesson." ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "First, we include the necessary headers for the ```MueLu::Hierarchy```" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_MutuallyExclusiveTime.hpp\"\n", "#include \"MueLu_TimeMonitor.hpp\"" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_Hierarchy.hpp\"" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then, we create a new ```MueLu::Hierarchy``` with standard multigrid parameters." ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::Hierarchy &) @0x7f4a182be950\n" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::Hierarchy H(xMat);" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "We set the maximum size for the coarsest level to $10$ and increase the verbosity output level to get more screen output." ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "H.setDefaultVerbLevel(Teuchos::VERB_HIGH);\n", "H.SetMaxCoarseSize(10); " ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "##### SmootherFactory" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "Then, we create a ```MueLu::SmootherFactory``` which serves as general interface for the smoothers and direct solvers for the multigrid hierarchy." ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_SmootherFactory.hpp\"" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::SmootherFactory &) @0x7f4a182bead8\n" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::SmootherFactory smootherFact;\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "We fill the *SmootherFactory* with a concrete smoother. In this case, we use the ```trSmoo``` smoother object for pre- and post smoothing. One could also provide two distinct smoothers/solvers for pre-/postsmoothing or use the ```MueLu::NoSmoother```" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::TrilinosSmoother &) @0x7f4a182bec50\n" ] }, "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::TrilinosSmoother trSmoo(\"RELAXATION\", subSmooParams);" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smootherFact.SetSmootherPrototypes(Teuchos::rcpFromRef(trSmoo));" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "#### The FactoryManager" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "Then, we create a *FactoryManager* which manages the default instances for each variable.\n", "More details on the *FactoryManager* can be found in the *Hierarchy* lesson." ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#include \"MueLu_FactoryManager.hpp\"" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(MueLu::FactoryManager &) @0x7f4a182bee30\n" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MueLu::FactoryManager M;" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "deletable": true, "editable": true }, "source": [ "Use the ```SetFactory``` routine to register the ```SmootherFactory``` to produce the ```Smoother``` variable (and the ```DirectSolver```).\n", "
Note:\n", "For real examples one would use a different factory providing an actual direct solver from the *Amesos* package, of course.\n", "
" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false, "deletable": true, "editable": true }, "outputs": [ { "data": { "text/plain": [ "(void) @0x7ffda227d1d0\n" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "M.SetFactory(\"Smoother\", Teuchos::rcpFromRef(smootherFact));\n", "M.SetFactory(\"CoarseSolver\", Teuchos::rcpFromRef(smootherFact));\n" ] }, { "cell_type": "markdown", "metadata": { "deletable": true, "editable": true }, "source": [ "With these factories set in the ```FactoryManager``` you can create a hierarchy (assuming that all other necessary variables and Factories are properly defined; please see the corresponding notebook on ```MueLu::Hierarchy```)." ] } ], "metadata": { "kernelspec": { "display_name": "C++11", "language": "", "name": "cling-cpp11" }, "language_info": { "codemirror_mode": "text/x-c++src", "file_extension": ".c++", "mimetype": " text/x-c++src", "name": "c++" } }, "nbformat": 4, "nbformat_minor": 0 }