mixed: document select structs and functions

opm/simulators/linalg/mixed/README.md

@@ -1,17 +1,19 @@
 # Mixed-precision linear solvers
 This folder contains mixed-precision building blocks for Krylov subspace methods
-and a highly optimized mixed-precision implementation of ILU0 preconditioned bicgstab.
-Hopefully, this will inspire the exploration of mixed-precision algorithms in OPM.
+and a highly optimized mixed-precision implementation of ILU0 and DILU preconditioned
+bicgstab. Hopefully, this will inspire the exploration of mixed-precision algorithms
+in OPM.
 
 The initial implementation is specialized for 3x3 block-sparse matrices due to their
 importance in reservoir simulation and restricted to serial runs. Extending the work
 to parallel runs is a matter of extending OPM's parallel infrastructure and should be
 relatively straight-forward.
 
-The mixed-precision solver is selected by the command-line option `--linear-solver=mixed-ilu0`.
-The command-line option `--matrix-add-well-contributions=true` must also be set as the
-mixed-precision solver operates directly on block-sparse matrices, not on linear operators as
-other OPM solvers do. For convenience, a wrapper similar to the one below can be used.
+The mixed-precision solver is selected by the command-line options `--linear-solver=mixed-ilu0`
+or `--linear-solver=mixed-dilu`. The command-line option `--matrix-add-well-contributions=true`
+must also be set as the mixed-precision solver operates directly on block-sparse matrices, not
+on linear operators as other OPM solvers do. For convenience, a wrapper similar to the one below
+can be used.
 
 ``` bash
 OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \

opm/simulators/linalg/mixed/bslv.h

@@ -10,19 +10,28 @@ extern "C" {
 
 #include <stdbool.h>
 
-// Solver memory struct
+/*!
+ * @brief Linear solver memory.
+ */
 typedef
 struct bslv_memory
 {
+    // perform dilu (ilu0) of true (False)
     bool use_dilu;
 
+    // tolerance
     double tol;
+    // maximum number of iterations
     int    max_iter;
+    // convergencey history
     double *e;
 
+    // size of linear system
     int n;
+    // pointer to temporary arrays
     double **dtmp;
 
+    // pointer to preconditioner
     prec_t *P;
 }
 bslv_memory;

opm/simulators/linalg/mixed/bsr.h

@@ -4,18 +4,28 @@
 extern "C" {
 #endif
 
-// bsr matrix struct
+/*!
+ * @brief Mixed-precision bsr matrix.
+ */
 typedef
 struct bsr_matrix
 {
+    // number or rows
     int nrows;
+    // number or columns
     int ncols;
+    // number or non-zero blocks
     int nnz;
+    // block-size
     int b;
 
+    // pointer ot row offsets
     int *rowptr;
+    // pointer to column indices
     int *colidx;
+    // pointer to double-precision values
     double *dbl;
+    // pointer to single-precision values
     float *flt;
 
 } bsr_matrix;

opm/simulators/linalg/mixed/prec.c

@@ -122,6 +122,13 @@ void prec_init(prec_t *P, bsr_matrix const *A)
     P->noffsets=count;
 }
 
+/**
+ * @brief Matrix-matrix multiplication of 3x3 matrices.
+ *
+ * @param C Pointer to output matrix.
+ * @param A Pointer to left input matrix.
+ * @param B Pointer to right input matrix.
+ */
 void mat3_matmul(double *C, const double *A, const double *B)
 {
     // assume 3x3 column-major matrices
@@ -141,6 +148,13 @@ void mat3_matmul(double *C, const double *A, const double *B)
     for(int k=0;k<9;k++) C[k]=M[k];
 }
 
+/**
+ * @brief Fused multiply-subtract for 3x3 matrices.
+ *
+ * @param C Pointer to output matrix.
+ * @param A Pointer to left input matrix.
+ * @param B Pointer to right input matrix.
+ */
 void mat3_matfms(double *C, const double *A, const double *B)
 {
     // assume 3x3 column-major matrices
@@ -160,7 +174,12 @@ void mat3_matfms(double *C, const double *A, const double *B)
     for(int k=0;k<9;k++) C[k]-=M[k];
 }
 
-
+/**
+ * @brief Matrix invese for 3x3 matrix.
+ *
+ * @param invA Pointer to inverse matrix.
+ * @param    A Pointer to input matrix.
+ */
 void mat3_inv(double *invA, const double *A)
 {
     double M[9];
@@ -178,11 +197,25 @@ void mat3_inv(double *invA, const double *A)
     for(int k=0;k<9;k++) invA[k]=M[k]/detA;
 }
 
+/**
+ * @brief vector copy of 9-element vectors.
+ *
+ * @param y Pointer to output vector.
+ * @param x Pointer to input vector.
+ */
 inline void vec_copy9(double *y, double const *x)
 {
     for(int i=0;i<9;i++) y[i]=x[i];
 }
 
+/**
+ * @brief In-place right matrix-matrix multiplication for 3x3 matrices.
+ *
+ * @note Function is optimized with avx2 intrinsics
+ *
+ * @param A Pointer to left input and output matrix.
+ * @param B Pointer to right input matrix.
+ */
 void mat3_rmul(double *A, double const *B)
 {
     // load left hand matrix
@@ -211,6 +244,14 @@ void mat3_rmul(double *A, double const *B)
     }
 }
 
+/**
+ * @brief In-place left matrix-matrix multiplication for 3x3 matrices.
+ *
+ * @note Function is optimized with avx2 intrinsics
+ *
+ * @param A Pointer to left input  matrix.
+ * @param B Pointer to right input and output matrix.
+ */
 void mat3_lmul(double const *A, double *B)
 {
     // load left hand matrix
@@ -239,6 +280,15 @@ void mat3_lmul(double const *A, double *B)
     }
 }
 
+/**
+ * @brief Fused multiply-subtract for 3x3 matrices.
+ *
+ * @note Function is optimized with avx2 intrinsics
+ *
+ * @param C Pointer to output matrix.
+ * @param A Pointer to left input matrix.
+ * @param B Pointer to right input matrix.
+ */
 void mat3_vfms(double *C, double const *A, double const *B)
 {
     // load left hand matrix
@@ -408,7 +458,6 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A)
             next=P->offsets[++idx][0];
         }
 
-
         for(int k=L->rowptr[i];k<L->rowptr[i+1];k++)
         {
             //scale row i of U
@@ -484,8 +533,6 @@ void prec_mapply3c(prec_t *restrict P, double *x)
             __m256d vxj = _mm256_loadu_pd(xj);
             __m256d vz = (vxj - vA[0]) - (vA[1] + vA[2]);
 
-            //vz =_mm256_blend_pd(vxj,vz,0x7);  // 4th element unchanged
-            //_mm256_storeu_pd(xj,vz);
             double z[4];
             _mm256_store_pd(z,vz);
             for(int n=0;n<3;n++) xj[n]=z[n];
@@ -498,8 +545,6 @@ void prec_mapply3c(prec_t *restrict P, double *x)
         vA[2] = _mm256_cvtps_pd(_mm_loadu_ps(A+6))*vx[2]; //0b01010101
         __m256d vz = vA[0] + vA[1] + vA[2];
 
-        //vz =_mm256_blend_pd(vxi,vz,0x7);  // 4th element unchanged
-        //_mm256_storeu_pd(xi,vz);
         double z[4];
         _mm256_store_pd(z,vz);
         for(int k=0;k<3;k++) xi[k]=z[k];
@@ -524,10 +569,6 @@ void prec_mapply3c(prec_t *restrict P, double *x)
         __m256d vz = (vxi - vA[0]) - (vA[1] + vA[2]);
         vz =_mm256_blend_pd(vxi,vz,0x7);  // 4th element unchanged
         _mm256_storeu_pd(xi,vz);
-
-        //double z[4];
-        //_mm256_store_pd(z,vz);
-        //for(int k=0;k<3;k++) xi[k]=z[k];
     }
 }
 
@@ -565,8 +606,6 @@ void prec_dapply3c(prec_t *restrict P, double *x)
             __m256d vxj = _mm256_loadu_pd(xj);
             __m256d vz = (vxj - vA[0]) - (vA[1] + vA[2]);
 
-            //vz =_mm256_blend_pd(vxj,vz,0x7);  // 4th element unchanged
-            //_mm256_storeu_pd(xj,vz);
             double z[4];
             _mm256_store_pd(z,vz);
             for(int kk=0;kk<3;kk++) xj[kk]=z[kk];
@@ -579,8 +618,6 @@ void prec_dapply3c(prec_t *restrict P, double *x)
         vA[2] = _mm256_loadu_pd(A+6)*vx[2]; //0b01010101
         __m256d vz = vA[0] + vA[1] + vA[2];
 
-        //vz =_mm256_blend_pd(vxi,vz,0x7);  // 4th element unchanged
-        //_mm256_storeu_pd(xi,vz);
         double z[4];
         _mm256_store_pd(z,vz);
         for(int k=0;k<3;k++) xi[k]=z[k];
@@ -605,10 +642,6 @@ void prec_dapply3c(prec_t *restrict P, double *x)
         __m256d vz = (vxi - vA[0]) - (vA[1] + vA[2]);
         vz =_mm256_blend_pd(vxi,vz,0x7);  // 4th element unchanged
         _mm256_storeu_pd(xi,vz);
-
-        //double z[4];
-        //_mm256_store_pd(z,vz);
-        //for(int k=0;k<3;k++) xi[k]=z[k];
     }
 }
 

opm/simulators/linalg/mixed/prec.h

@@ -6,14 +6,21 @@ extern "C" {
 
 #include "bsr.h"
 
-// preconditioner struct
+/*!
+ * @brief Preconditioner struct.
+ */
 typedef
 struct prec_t
 {
+    // lower triangular factor
     bsr_matrix *L;
+    // diagonal factor
     bsr_matrix *D;
+    // upper triangular factor
     bsr_matrix *U;
+    // number of off-diagonal ilu0 updates
     int noffsets;
+    // triplets uniquely identifying off-diagonal ilu0 updates
     int(*offsets)[3];
 }
 prec_t;

opm/simulators/linalg/mixed/vec.h

@@ -5,11 +5,25 @@ extern "C" {
 #endif
 
 
+/**
+ * @brief Copy vector.
+ *
+ * @param y Pointer to output vector.
+ * @param x Pointer to inout vector.
+ * @param n Length of vectors.
+ */
 void vec_copy(double *y, double const * x, int n)
 {
     for(int i=0;i<n;i++) y[i]=x[i];
 }
 
+/**
+ * @brief Fill vector.
+ *
+ * @param y Pointer to output vector.
+ * @param x Fill value.
+ * @param n Length of vector.
+ */
 void vec_fill(double *y, double x, int n)
 {
     for(int i=0;i<n;i++) y[i]=x;