{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I recently discovered [cling](https://root.cern.ch/drupal/content/cling) after I saw this Tweet\n",
    "\n",
    "<blockquote class=\"twitter-tweet\" lang=\"en\"><p lang=\"en\" dir=\"ltr\">I made a C++ kernel for Jupyter for some reason: <a href=\"https://t.co/7ApuZyeR0r\">https://t.co/7ApuZyeR0r</a></p>&mdash; Min RK (@minrk) <a href=\"https://twitter.com/minrk/status/600762804161933312\">May 19, 2015</a></blockquote>\n",
    "<script async src=\"//platform.twitter.com/widgets.js\" charset=\"utf-8\"></script>\n",
    "\n",
    "From the cling website, \"Cling is an interactive C++ interpreter, built on the top of LLVM and Clang libraries. Its advantages over the standard interpreters are that it has command line prompt and uses just-in-time (JIT) compiler for compilation.\"  As I write (or mostly teach through dual-coding with students) quite a bit of C++ and have already been utilizing the [IPython](http://ipython.org/)/[Jupyter](https://jupyter.org/) notebook as a teaching tool for sometime, I was immediately interested in using [MinRK](https://github.com/minrk)'s [clingkernel](https://github.com/minrk/clingkernel) for this purpose as well as exploring and testing out my own ideas.  As I learned more about cling, I found most of the examples to be very trivial and those that inluded the calling/loading of external functions typically only utilized the functions from the standard library.  I was eager to give it try on something a little more complex, and because most of the code we write in my group is heavily dependent on [Trilinos](http://trilinos.org/), I thought I would attempt to solve a simple linear algebra system of the form\n",
    "\n",
    "$$\n",
    "A \\vec{x} = \\vec{b}\n",
    "$$\n",
    "\n",
    "using [AztecOO](http://trilinos.org/packages/epetra/) solvers with [Epetra](http://trilinos.org/packages/aztecoo/) datastructures from Trilinos (at least AztecOO and Epetra).  This is an adaptation of the example code provided in Figure 1 of the [AztecOO User's guide](http://trilinos.org/oldsite/packages/aztecoo/AztecOOUserGuide.pdf).\n",
    "\n",
    "Of course, to reproduce the code below, you will need to [install cling](https://github.com/vgvassilev/cling) and Trilinos.  Additionally, if you want to use MinRK's [clingkernel](https://github.com/minrk/clingkernel) you need to install the [Jupyter notebook](https://github.com/jupyter/notebook#jupyter-notebook) from the `master` branch, and install the kernel.  With that, you should be able to <a href=\"http://johntfoster.github.io/posts/interactive-c++-for-hpc.ipynb\">download this notebook</a> and execute it with few modifications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cling provides a few metaprocessing commands to perform actions that you would normal provide as arguments to the build system, such as providing locations for include files, loading shared libraries, etc.  The documentation for the metaprocessing commands are\n",
    "\n",
    "````\n",
    ".x test.cpp – load a file (if the file defines void test() it will be executed)\n",
    ".L libname – load a libary or file\n",
    ".x filename.cxx – loads filename and calls void filename() if defined\n",
    ".I path – adds an include path\n",
    ".printAST – shows the abstract syntax tree after each processed entity\n",
    ".q – exit cling\n",
    ".help – display a brief description\n",
    "````\n",
    "\n",
    "We'll start by loading the location of the installed header files in my Trilinos build (here I only built AztecOO and Epetra in serial to keep it simple).  Cling utilizes the Clang compiler which is pretty picky, so I had to edit the header files in a few locations where `#include` statments where provided with angle brackets `< >` for local includes, where the compiler wanted double quotes.  In other words, I had to change\n",
    "\n",
    "````\n",
    "#include <foo.h>\n",
    "````\n",
    "\n",
    "to \n",
    "\n",
    "````\n",
    "#include \"foo.h\"\n",
    "````\n",
    "\n",
    "in a few of the Trilinos headers.  So the metaprocessing command to tell the system where the header files are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    ".I /usr/local/trilinos-serial/include/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we provide the include statememts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "#include \"AztecOO.h\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "#include \"AztecOO_Version.h\"\n",
    "#include \"Epetra_SerialComm.h\"\n",
    "#include \"Epetra_Map.h\"\n",
    "#include \"Epetra_Vector.h\"\n",
    "#include \"Epetra_CrsMatrix.h\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the metaprocessing to load the shared libraries.  Since I'm on a Mac, I had to define the environment variable \n",
    "\n",
    "````bash\n",
    "export DYLD_LIBRARY_PATH=/usr/local/trilinos-serial/lib:$DYLD_LIBRARY_PATH\n",
    "````\n",
    "\n",
    "before launching the notebook/cling so the system knows where to find the libraries at runtime.  On a Linux machine, this variable would be `LD_LIBRARY_PATH`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    ".L libepetra\n",
    ".L libaztecOO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we start in with the actual Epetra/AztecOO code.  Since it's not apparent to me whether it's possible for `cling` to support [MPI](http://www.open-mpi.org/) processes, we'll just provide a serial implementation.  We start with some code that instantiates the Epetra serial communicator and provides some information about the version of Trilinos we are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "AztecOO in Trilinos 12.1 (Dev)\r\n",
      "\r\n",
      "\r\n",
      "Epetra::Comm::Processor 0 of 1 total processors.\r\n"
     ]
    }
   ],
   "source": [
    "Epetra_SerialComm Comm;\n",
    "  \n",
    "int NumMyElements = 1000;\n",
    "\n",
    "if (Comm.MyPID()==0)\n",
    "    std::cout << AztecOO_Version() << std::endl << std::endl;\n",
    "\n",
    "  std::cout << Comm <<std::endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice the immediate output to screen.  Remember this is C++!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we construct a Map that puts same number of equations on each processor.  Of course, we only have one processor here, but the code is generic and could also be run in parallel if multiple processors were available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "Epetra_Map Map(-1, NumMyElements, 0, Comm);\n",
    "int NumGlobalElements = Map.NumGlobalElements();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we instantiate the column-row sparse matrix $A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "Epetra_CrsMatrix A(Copy, Map, 3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we fill $A$ to be a finite-difference like operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "?   \r\n"
     ]
    }
   ],
   "source": [
    "double negOne = -1.0;\n",
    "double posTwo = 2.0;\n",
    "for (int i=0; i<NumMyElements; i++) {\n",
    "    int GlobalRow = A.GRID(i); int RowLess1 = GlobalRow - 1; int RowPlus1 = GlobalRow + 1;\n",
    "\n",
    "    if (RowLess1!=-1) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowLess1);\n",
    "    if (RowPlus1!=NumGlobalElements) A.InsertGlobalValues(GlobalRow, 1, &negOne, &RowPlus1);\n",
    "    A.InsertGlobalValues(GlobalRow, 1, &posTwo, &GlobalRow);\n",
    "};"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we call the `FillComplete()` method to optimize storage of the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": []
    }
   ],
   "source": [
    "A.FillComplete();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we instantiate $\\vec{x}$ and the right-hand-side vector $\\vec{b}$.  Also, we fill $\\vec{b}$ with random numbers for our solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n"
     ]
    }
   ],
   "source": [
    "Epetra_Vector x(Map);\n",
    "Epetra_Vector b(Map);\n",
    "\n",
    "b.Random();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we instantiate the `Epetra_Problem` and the Aztec solver from that instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n"
     ]
    }
   ],
   "source": [
    "Epetra_LinearProblem problem(&A, &x, &b);\n",
    "\n",
    "AztecOO solver(problem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next is the actual solution step.  A GMRES solver is used for 10 iterations.  This won't be enough for convergence of the solver, but this is just for demonstration purposes and I would like to keep the output short."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\r\n",
      "\t\t*******************************************************\r\n",
      "\t\t***** Problem: Epetra::CrsMatrix\r\n",
      "\t\t***** Preconditioned GMRES solution\r\n",
      "\t\t***** 1 step block Jacobi\r\n",
      "\t\t***** No scaling\r\n",
      "\t\t*******************************************************\r\n",
      "\r\n",
      "                iter:    0           residual = 1.000000e+00\r\n",
      "                iter:    1           residual = 5.832906e-01\r\n",
      "                iter:    2           residual = 4.547092e-01\r\n",
      "                iter:    3           residual = 3.831399e-01\r\n",
      "                iter:    4           residual = 3.366700e-01\r\n",
      "                iter:    5           residual = 3.043737e-01\r\n",
      "                iter:    6           residual = 2.789514e-01\r\n",
      "                iter:    7           residual = 2.592427e-01\r\n",
      "                iter:    8           residual = 2.449521e-01\r\n",
      "                iter:    9           residual = 2.328859e-01\r\n",
      "                iter:   10           residual = 2.223153e-01\r\n",
      "\r\n",
      "\r\n",
      "\t***************************************************************\r\n",
      "\r\n",
      "\tWarning: maximum number of iterations exceeded\r\n",
      "\twithout convergence\r\n",
      "\r\n",
      "\tSolver:\t\t\tgmres\r\n",
      "\tnumber of iterations:\t10\r\n",
      "\r\n",
      "\tRecursive residual =  3.9447e+00\r\n",
      "\r\n",
      "\tCalculated Norms\t\t\t\tRequested Norm\r\n",
      "\t--------------------------------------------\t--------------\r\n",
      "\r\n",
      "\t||r||_2 / ||r0||_2:\t\t2.223153e-01\t1.000000e-02\r\n",
      "\r\n",
      "\t***************************************************************\r\n",
      "\r\n",
      "\r\n",
      "\r\n",
      "\t\tSolution time: 0.000680 (sec.)\r\n",
      "\t\ttotal iterations: 10\r\n",
      "\r\n",
      "Solver performed 10 iterations.\r\n",
      "Norm of true residual = 3.94466\r\n"
     ]
    }
   ],
   "source": [
    "solver.SetAztecOption(AZ_precond, AZ_Jacobi);\n",
    "solver.Iterate(10, 1.0E-2);\n",
    "\n",
    "std::cout << \"Solver performed \" << solver.NumIters() << \" iterations.\" << std::endl;\n",
    "std::cout << \"Norm of true residual = \" << solver.TrueResidual() << std::endl;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To demonstrate that is worked, we'll print out the first 10 values of the solution vector $\\vec{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.288916\r\n",
      "-0.723977\r\n",
      "-0.641509\r\n",
      "0.309992\r\n",
      "0.815517\r\n",
      "0.904473\r\n",
      "0.0381509\r\n",
      "-0.675333\r\n",
      "-1.5226\r\n",
      "-3.14536\r\n"
     ]
    }
   ],
   "source": [
    "for (int i = 0; i < 10; i++){\n",
    "    std::cout << x[i] << std::endl; \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's pretty amazing to be able to do interactive C++ developement.  If we could get this to work in parallel, it would be a leap forward in the ability to do introspective high-performance scientific computing."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "C++",
   "language": "",
   "name": "cling"
  },
  "language_info": {
   "codemirror_mode": "clike",
   "file_extension": ".c++",
   "mimetype": " text/x-c++src",
   "name": "c++"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}