/*

* basiclu_factorize.c

*

* Copyright (C) 2016-2018 ERGO-Code

*

*/

#include "lu_internal.h"

#include "lu_timer.h"

lu_int basiclu_factorize

(

lu_int istore[],

double xstore[],

lu_int Li[],

double Lx[],

lu_int Ui[],

double Ux[],

lu_int Wi[],

double Wx[],

const lu_int Bbegin[],

const lu_int Bend[],

const lu_int Bi[],

const double Bx[],

lu_int c0ntinue

)

{

struct lu this;

lu_int status;

double tic[2], elapsed, factor_cost;

lu_tic(tic);

status = lu_load(&this, istore, xstore, Li, Lx, Ui, Ux, Wi, Wx);

if (status != BASICLU_OK)

return status;

if (! (Li && Lx && Ui && Ux && Wi && Wx && Bbegin && Bend && Bi && Bx))

{

status = BASICLU_ERROR_argument_missing;

return lu_save(&this, istore, xstore, status);

}

if (!c0ntinue)

{

lu_reset(&this);

this.task = SINGLETONS;

}

/* continue factorization */

switch (this.task)

{

case SINGLETONS: goto singletons;

case SETUP_BUMP: goto setup_bump;

case FACTORIZE_BUMP: goto factorize_bump;

case BUILD_FACTORS: goto build_factors;

}

status = BASICLU_ERROR_invalid_call;

return lu_save(&this, istore, xstore, status);

/*

* Each of the following four parts of the factorization calls a routine

* lu_do_something() which may request reallocation. In this case return

* to the caller immediately, keeping the entry point in this.task.

*/

singletons:

this.task = SINGLETONS;

status = lu_singletons(&this, Bbegin, Bend, Bi, Bx);

if (status != BASICLU_OK)

goto return_to_caller;

setup_bump:

this.task = SETUP_BUMP;

status = lu_setup_bump(&this, Bbegin, Bend, Bi, Bx);

if (status != BASICLU_OK)

goto return_to_caller;

factorize_bump:

this.task = FACTORIZE_BUMP;

status = lu_factorize_bump(&this);

if (status != BASICLU_OK)

goto return_to_caller;

build_factors:

this.task = BUILD_FACTORS;

status = lu_build_factors(&this);

if (status != BASICLU_OK)

goto return_to_caller;

/* factorization successfully finished */

this.task = NO_TASK;

this.nupdate = 0; /* make factorization valid */

this.ftran_for_update = -1;

this.btran_for_update = -1;

this.nfactorize++;

this.condestL = lu_condest(this.m, this.Lbegin, this.Lindex, this.Lvalue,

NULL, this.p, 0, this.work1,

&this.normL, &this.normestLinv);

this.condestU = lu_condest(this.m, this.Ubegin, this.Uindex, this.Uvalue,

this.row_pivot, this.p, 1, this.work1,

&this.normU, &this.normestUinv);

/* measure numerical stability of the factorization */

lu_residual_test(&this, Bbegin, Bend, Bi, Bx);

/*

* factor_cost is a deterministic measure of the factorization cost.

* The parameters have been adjusted such that (on my computer)

* 1e-6 * factor_cost =~ time_factorize.

*

* update_cost measures the accumulated cost of updates/solves compared

* to the last factorization. It is computed from

*

* update_cost = update_cost_numer / update_cost_denom.

*

* update_cost_denom is fixed here.

* update_cost_numer is zero here and increased by solves/updates.

*

*/

factor_cost =

0.04 * this.m +

0.07 * this.matrix_nz +

0.20 * this.bump_nz +

0.20 * this.nsearch_pivot +

0.008 * this.factor_flops;

this.update_cost_denom = factor_cost * 250;

#if 0

elapsed = this.time_factorize + lu_toc(tic);

printf(" 1e-6 * factor_cost / time_factorize: %.2f\n",

1e-6 * factor_cost / elapsed);

#endif

if (this.rank < this.m)

status = BASICLU_WARNING_singular_matrix;

return_to_caller:

elapsed = lu_toc(tic);

this.time_factorize += elapsed;

this.time_factorize_total += elapsed;

return lu_save(&this, istore, xstore, status);

}