From ea230f63121cec85a81c4796a4e3ce88979860a3 Mon Sep 17 00:00:00 2001 From: Aznaveh Date: Wed, 19 Jun 2024 20:14:28 -0500 Subject: [PATCH] adding some stats --- ParU/Demo/paru_simple.cpp | 4 ++++ ParU/Source/ParU_Factorize.cpp | 20 +++++++++++++++++++- ParU/Source/ParU_Get.cpp | 6 +++--- ParU/Source/paru_internal.hpp | 3 +++ 4 files changed, 29 insertions(+), 4 deletions(-) diff --git a/ParU/Demo/paru_simple.cpp b/ParU/Demo/paru_simple.cpp index f27895eb7..ba63f6bb1 100644 --- 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 index e194b5f0b..2297d40d0 100644 --- 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 index d82716f87..4dde9c688 100644 --- 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 index f0445aa9c..59c028b43 100644 --- 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 } ;