adding some stats

DrTimothyAldenDavis · Jun 20, 2024 · ea230f6 · ea230f6
1 parent 8360819
commit ea230f6
Show file tree

Hide file tree

Showing 4 changed files with 29 additions and 4 deletions.
diff --git a/ParU/Demo/paru_simple.cpp b/ParU/Demo/paru_simple.cpp
@@ -70,7 +70,11 @@ int main(int argc, char **argv)
     std::cout << "Input matrix is " << n << "x" << n <<
         " nnz = " << anz << std::endl;
     OK (ParU_Factorize(A, Sym, &Num, Control), "numeric factorization") ;
+    int64_t unz, lnz ;
+    OK (ParU_Get (Sym, Num, PARU_GET_LNZ, &lnz, Control), "lnz") ;
+    OK (ParU_Get (Sym, Num, PARU_GET_UNZ, &unz, Control), "unz") ;
     std::cout << "ParU: factorization was successful." << std::endl;
+    std::cout << "nnz(L) = " << lnz << ", nnz(U) = " << unz << std::endl;
 
     //~~~~~~~~~~~~~~~~~~~ Computing the residual, norm(b-Ax) ~~~~~~~~~~~~~~~~~~~
     b = (double *)malloc(n * sizeof(double));

diff --git a/ParU/Source/ParU_Factorize.cpp b/ParU/Source/ParU_Factorize.cpp
@@ -282,6 +282,9 @@ ParU_Info ParU_Factorize
 #endif
     int64_t max_rc = 0, max_cc = 0;
     double min_udiag = 1, max_udiag = -1;  // not to fail for nf ==0
+    int64_t nnzL=0;
+    int64_t nnzU=0;
+    double sfc = 0.0; //simple flop count
 
     // using the first value of the first front just to initialize
     if (nf > 0)
@@ -307,6 +310,12 @@ ParU_Info ParU_Factorize
                 max_rc = std::max(max_rc, rowCount);
                 max_cc = std::max(max_cc, colCount + fp);
                 double *X = LUs[f].p;
+                nnzL += fp*rowCount - fp*(fp-1)/2;
+                nnzU += fp*colCount + fp*(fp+1)/2;
+                sfc+= (4*fp*fp*fp-3*fp*fp-fp)/6;
+                sfc+= (fp*(fp-1)/2+fp)*(rowCount-fp);
+                sfc+= (fp*(fp-1)/2)*(colCount);
+                sfc+= (rowCount-fp)*(colCount);
                 for (int64_t i = 0; i < fp; i++)
                 {
                     double udiag = fabs(X[rowCount * i + i]);
@@ -332,13 +341,19 @@ ParU_Info ParU_Factorize
                 max_rc = std::max(max_rc, rowCount);
                 max_cc = std::max(max_cc, colCount + fp);
             }
-
             for (int64_t f = 0; f < nf; f++)
             {
                 int64_t rowCount = Num->frowCount[f];
+                int64_t colCount = Num->fcolCount[f];
                 int64_t col1 = Super[f];
                 int64_t col2 = Super[f + 1];
                 int64_t fp = col2 - col1;
+                nnzL += fp*rowCount - fp*(fp-1)/2;
+                nnzU += fp*colCount + fp*(fp+1)/2;
+                sfc+= (4*fp*fp*fp-3*fp*fp-fp)/6;
+                sfc+= (fp*(fp-1)/2+fp)*(rowCount-fp);
+                sfc+= (fp*(fp-1)/2)*(colCount);
+                sfc+= (rowCount-fp)*(colCount);
                 double *X = LUs[f].p;
                 #pragma omp parallel for reduction(min:min_udiag)   \
                     reduction(max: max_udiag)                       \
@@ -361,6 +376,9 @@ ParU_Info ParU_Factorize
     Num->min_udiag = min_udiag;
     Num->max_udiag = max_udiag;
     Num->rcond = min_udiag / max_udiag;
+    Num->nnzL = nnzL;
+    Num->nnzU = nnzU;
+    Num->sfc= sfc;
 #ifndef NTIME
     double time = PARU_OPENMP_GET_WTIME;
     time -= my_start_time;

diff --git a/ParU/Source/ParU_Get.cpp b/ParU/Source/ParU_Get.cpp
@@ -81,12 +81,12 @@ ParU_Info ParU_Get              // get an int64_t scalar or array from Sym/Num
 
         case PARU_GET_LNZ:
             if (!Num || Num->sym_m != n) return (PARU_INVALID) ;
-            (*result) = 0 ;     // FIXME NOW: get nnz(L)
+            (*result) = (int64_t) Num->nnzL ;
             break ;
 
         case PARU_GET_UNZ:
             if (!Num || Num->sym_m != n) return (PARU_INVALID) ;
-            (*result) = 0 ;     // FIXME NOW: get nnz(U)
+            (*result) = (int64_t) Num->nnzU ;
             break ;
 
         case PARU_GET_P:
@@ -150,7 +150,7 @@ ParU_Info ParU_Get              // get a double scalar or array from Sym/Num
     switch (field)
     {
         case PARU_GET_FLOP_COUNT:
-            (*result) = 0 ;     // FIXME NOW: get flop count
+            (*result) = (double) Num->sfc ;     
             break ;
 
         case PARU_GET_RCOND_ESTIMATE:

diff --git a/ParU/Source/paru_internal.hpp b/ParU/Source/paru_internal.hpp
@@ -241,6 +241,9 @@ struct ParU_Numeric_struct
     double rcond;       // rough estimate of the reciprocal condition number
     double min_udiag;   // min (abs (diag (U)))
     double max_udiag;   // max (abs (diag (U)))
+    int64_t nnzL;       //nnz of L
+    int64_t nnzU;       //nnz of U
+    double sfc; //simple flop count
     ParU_Info res;  // returning value of numeric phase
 } ;