Skip to content

Commit

Permalink
.
Browse files Browse the repository at this point in the history
  • Loading branch information
DrTimothyAldenDavis committed Feb 20, 2024
1 parent 3b8a459 commit fa147f6
Show file tree
Hide file tree
Showing 7 changed files with 4 additions and 110 deletions.
10 changes: 0 additions & 10 deletions SPEX/SPEX_Cholesky/Source/SPEX_cholesky_backslash.c
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,6 @@ SPEX_info SPEX_cholesky_backslash
const SPEX_options option // Command options (Default if NULL)
)
{
HERE

SPEX_info info;
// SPEX must be initialized
Expand Down Expand Up @@ -112,11 +111,9 @@ HERE

bool is_symmetric ;
SPEX_CHECK( SPEX_determine_symmetry(&is_symmetric, A, option) );
HERE
if (!is_symmetric)
{
SPEX_FREE_WORKSPACE ;
HERE
return SPEX_NOTSPD ;
}

Expand All @@ -125,22 +122,19 @@ HERE
//--------------------------------------------------------------------------

SPEX_CHECK( spex_cholesky_preorder(&S, A, option) );
HERE

//--------------------------------------------------------------------------
// Permute matrix A, that is apply the row/column ordering from the
// symbolic analysis step to get the permuted matrix PAP.
//--------------------------------------------------------------------------

SPEX_CHECK( spex_cholesky_permute_A(&PAP, A, true, S) );
HERE

//--------------------------------------------------------------------------
// Symbolic Analysis: compute the elimination tree of PAP
//--------------------------------------------------------------------------

SPEX_CHECK( spex_cholesky_symbolic_analysis(S, PAP, option) );
HERE

//--------------------------------------------------------------------------
// Factorization: Perform the REF Cholesky factorization of PAP.
Expand All @@ -149,7 +143,6 @@ HERE
//--------------------------------------------------------------------------

SPEX_CHECK( spex_cholesky_factor(&F, S, PAP, option) );
HERE

//--------------------------------------------------------------------------
// Solve: Solve Ax = b using the REF Cholesky factorization. That is,
Expand All @@ -159,7 +152,6 @@ HERE
//--------------------------------------------------------------------------

SPEX_CHECK( SPEX_cholesky_solve(&x, F, b, option) );
HERE

//--------------------------------------------------------------------------
// At this point x is stored as mpq_t. If the user desires the output
Expand All @@ -178,13 +170,11 @@ HERE
(*x_handle) = x2;
SPEX_matrix_free (&x, NULL);
}
HERE

//--------------------------------------------------------------------------
// Free all workspace and return success
//--------------------------------------------------------------------------

SPEX_FREE_WORKSPACE;
HERE
return SPEX_OK;
}
11 changes: 0 additions & 11 deletions SPEX/SPEX_Cholesky/Source/spex_cholesky_factor.c
Original file line number Diff line number Diff line change
Expand Up @@ -71,15 +71,13 @@ SPEX_info spex_cholesky_factor
// reminder of the appropriate data types
ASSERT(A->type == SPEX_MPZ);
ASSERT(A->kind == SPEX_CSC);
HERE

// Number of nonzeros in A
int64_t anz;
SPEX_CHECK( SPEX_matrix_nnz(&anz, A, option) );
ASSERT(anz > 0);

(*F_handle) = NULL ;
HERE

//--------------------------------------------------------------------------
// Declare and initialize workspace
Expand All @@ -90,7 +88,6 @@ HERE

// Allocate memory for the factorization
F = (SPEX_factorization) SPEX_calloc(1, sizeof(SPEX_factorization_struct));
HERE
if (F == NULL) return SPEX_OUT_OF_MEMORY;

// set factorization kind
Expand All @@ -99,20 +96,17 @@ HERE
// Allocate and set scale_for_A
SPEX_MPQ_INIT(F->scale_for_A);
SPEX_MPQ_SET(F->scale_for_A, A->scale);
HERE

// Inverse pivot ordering
F->Pinv_perm = (int64_t*) SPEX_malloc ( n*sizeof(int64_t) );
// row/column permutation, to be copied from S->P_perm
F->P_perm = (int64_t*) SPEX_malloc ( n*sizeof(int64_t) );
HERE
if (!(F->Pinv_perm) || !(F->P_perm))
{
// out of memory: free everything and return
SPEX_FREE_ALL;
return SPEX_OUT_OF_MEMORY;
}
HERE

// Copy row/column permutation from symbolic analysis to factorization
memcpy(F->P_perm, S->P_perm, n*sizeof(int64_t));
Expand All @@ -122,7 +116,6 @@ HERE
// factorization: up-looking or left-looking
//--------------------------------------------------------------------------

HERE
SPEX_factorization_algorithm algo = SPEX_OPTION_ALGORITHM(option);
switch(algo)
{
Expand All @@ -135,10 +128,8 @@ HERE
HERE
break;
case SPEX_CHOL_LEFT:
HERE
SPEX_CHECK( spex_cholesky_left_factor(&(F->L), &(F->rhos), S, A,
option) );
HERE
break;
default:
SPEX_FREE_ALL;
Expand All @@ -149,8 +140,6 @@ HERE
// Set outputs, return ok
//--------------------------------------------------------------------------

HERE
(*F_handle) = F ;
HERE
return SPEX_OK;
}
9 changes: 0 additions & 9 deletions SPEX/SPEX_Cholesky/Source/spex_cholesky_symbolic_analysis.c
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,6 @@ SPEX_info spex_cholesky_symbolic_analysis
{

SPEX_info info;
HERE

//--------------------------------------------------------------------------
// Check inputs
Expand All @@ -70,33 +69,25 @@ HERE
int64_t *c = NULL;

// Obtain elimination tree of A
HERE
SPEX_CHECK( spex_cholesky_etree(&S->parent, A) );

// Postorder the elimination tree of A
HERE
SPEX_CHECK( spex_cholesky_post(&post, S->parent, n) );

// Get the column counts of A
SPEX_CHECK( spex_cholesky_counts(&c, A, S->parent, post) );
HERE

// Set the column pointers of L
S->cp = (int64_t*) SPEX_malloc( (n+1)*sizeof(int64_t*));
HERE
if (S->cp == NULL)
{
SPEX_FREE_ALL;
HERE
return SPEX_OUT_OF_MEMORY;
}
HERE
SPEX_CHECK( spex_cumsum(S->cp, c, n));
HERE

// Set the exact number of nonzeros in L
S->lnz = S->cp[n];
SPEX_FREE_WORKSPACE;
HERE
return SPEX_OK;
}
26 changes: 0 additions & 26 deletions SPEX/SPEX_Cholesky/Source/spex_cholesky_up_factor.c
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,6 @@ SPEX_info spex_cholesky_up_factor
const SPEX_options option // command options
)
{
HERE

//--------------------------------------------------------------------------
// check inputs
Expand All @@ -89,7 +88,6 @@ HERE
ASSERT (rhos_handle != NULL);
(*L_handle) = NULL ;
(*rhos_handle) = NULL ;
HERE

//--------------------------------------------------------------------------
// Declare and initialize workspace
Expand All @@ -109,20 +107,17 @@ HERE
size_t size;

c = (int64_t*) SPEX_malloc(n*sizeof(int64_t));
HERE

// h is the history vector utilized for the sparse REF
// triangular solve algorithm. h serves as a global
// vector which is repeatedly passed into the triangular
// solve algorithm
h = (int64_t*) SPEX_malloc(n*sizeof(int64_t));
HERE

// xi serves as a global nonzero pattern vector. It stores
// the pattern of nonzeros of the kth column of L
// for the triangular solve.
xi = (int64_t*) SPEX_malloc(2*n*sizeof(int64_t));
HERE

if (!h || !xi || !c)
{
Expand All @@ -135,7 +130,6 @@ HERE
{
h[i] = -1;
}
HERE

//--------------------------------------------------------------------------
// Allocate and initialize the workspace x
Expand Down Expand Up @@ -165,14 +159,12 @@ HERE
// default size)
SPEX_CHECK (SPEX_matrix_allocate(&x, SPEX_DENSE, SPEX_MPZ, n, 1, n,
false, /* do not initialize the entries of x: */ false, option));
HERE

// Create rhos, a "global" dense mpz_t matrix of dimension n*1.
// As indicated with the second boolean parameter true, the mpz entries in
// rhos are initialized to the default size (unlike x).
SPEX_CHECK (SPEX_matrix_allocate(&(rhos), SPEX_DENSE, SPEX_MPZ, n, 1, n,
false, true, option));
HERE
printf ("estimate: % " PRId64"\n", estimate) ;
fflush (stdout) ;

Expand All @@ -182,7 +174,6 @@ HERE
// Allocate memory for entries of x to be estimate bits
SPEX_MPZ_INIT2(x->x.mpz[i], estimate);
}
HERE

//--------------------------------------------------------------------------
// Declare memory for L
Expand All @@ -196,15 +187,12 @@ HERE

SPEX_CHECK(SPEX_matrix_allocate(&(L), SPEX_CSC, SPEX_MPZ, n, n, S->lnz,
false, false, option));
HERE

// Set the column pointers of L
for (k = 0; k < n; k++)
{
L->p[k] = c[k] = S->cp[k];
}
HERE


//--------------------------------------------------------------------------
// Perform the up-looking factorization
Expand All @@ -214,7 +202,6 @@ HERE
// Iterations 0:n-1 (1:n in standard)
//--------------------------------------------------------------------------
SPEX_MPZ_SGN(&prev_sgn, x->x.mpz[0]);
HERE

for (k = 0; k < n; k++)
{
Expand All @@ -227,27 +214,23 @@ HERE
// If x[k] is nonzero choose it as pivot. Otherwise, the matrix is
// not SPD (indeed, it may even be singular).
SPEX_MPZ_SGN(&sgn, x->x.mpz[k]);
HERE
if (sgn != 0)
{
SPEX_MPZ_SET(rhos->x.mpz[k], x->x.mpz[k]);
}
else
{
// A is not symmetric positive definite
HERE
SPEX_FREE_ALL;
return SPEX_NOTSPD;
}
HERE

//----------------------------------------------------------------------
// Add the nonzeros (i.e. x) to L
//----------------------------------------------------------------------
int64_t p = 0;
for (j = top; j < n; j++)
{
HERE
// Obtain the row index of x[j]
jnew = xi[j];
if (jnew == k) continue;
Expand All @@ -258,34 +241,25 @@ HERE
// Place the i index of this nonzero. Should always be k because at
// iteration k, the up-looking algorithm computes row k of L
L->i[p] = k;
HERE

// Find the number of bits of x[j]
size = mpz_sizeinbase(x->x.mpz[jnew],2);
HERE

// GMP manual: Allocated size should be size+2
SPEX_MPZ_INIT2(L->x.mpz[p], size+2);
HERE

// Place the x value of this nonzero
SPEX_MPZ_SET(L->x.mpz[p],x->x.mpz[jnew]);
HERE
}
// Now, place L(k,k)
p = c[k]++;
L->i[p] = k;
HERE
size = mpz_sizeinbase(x->x.mpz[k], 2);
HERE
SPEX_MPZ_INIT2(L->x.mpz[p], size+2);
HERE
SPEX_MPZ_SET(L->x.mpz[p], x->x.mpz[k]);
HERE
}
// Finalize L->p
L->p[n] = S->lnz;
HERE

//--------------------------------------------------------------------------
// Free memory and set output
Expand Down
Loading

0 comments on commit fa147f6

Please sign in to comment.