1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
/* ========================================================================== */
/* === luflop mexFunction ================================================== */
/* ========================================================================== */
/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */
/* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack */
/* -------------------------------------------------------------------------- */
/*
f = luflop (L, U) ;
Given L and U, compute:
Lnz = full (sum (spones (L))) - 1 ;
Unz = full (sum (spones (U')))' - 1 ;
f = 2*Lnz*Unz + sum (Lnz) ;
without allocating O (lunz) space.
*/
#include "mex.h"
#include "matrix.h"
#ifndef TRUE
#define TRUE (1)
#endif
#ifndef FALSE
#define FALSE (0)
#endif
void mexFunction
(
int nlhs, /* number of left-hand sides */
mxArray *plhs [ ], /* left-hand side matrices */
int nrhs, /* number of right--hand sides */
const mxArray *prhs [ ] /* right-hand side matrices */
)
{
/* ---------------------------------------------------------------------- */
/* local variables */
/* ---------------------------------------------------------------------- */
double flop_count ;
double *pflop ;
int *Lp, *Li, *Up, *Ui, *Unz, n, k, row, col, p, Lnz_k, Unz_k ;
mxArray *Lmatrix, *Umatrix ;
/* ---------------------------------------------------------------------- */
/* get inputs L, U */
/* ---------------------------------------------------------------------- */
if (nrhs != 2)
{
mexErrMsgTxt ("Usage: f = luflop (L, U)") ;
}
Lmatrix = (mxArray *) prhs [0] ;
Umatrix = (mxArray *) prhs [1] ;
n = mxGetM (Lmatrix) ;
if (n != mxGetN (Lmatrix) || n != mxGetM (Umatrix) || n != mxGetN (Umatrix))
{
mexErrMsgTxt ("Usage: f = luflop (L, U) ; L and U must be square") ;
}
if (!mxIsSparse (Lmatrix) || !mxIsSparse (Umatrix))
{
mexErrMsgTxt ("Usage: f = luflop (L, U) ; L and U must be sparse") ;
}
Lp = mxGetJc (Lmatrix) ;
Li = mxGetIr (Lmatrix) ;
Up = mxGetJc (Umatrix) ;
Ui = mxGetIr (Umatrix) ;
Unz = (int *) mxMalloc (n * sizeof (int)) ;
/* ---------------------------------------------------------------------- */
/* count the nonzeros in each row of U */
/* ---------------------------------------------------------------------- */
for (row = 0 ; row < n ; row++)
{
Unz [row] = 0 ;
}
for (col = 0 ; col < n ; col++)
{
for (p = Up [col] ; p < Up [col+1] ; p++)
{
row = Ui [p] ;
Unz [row]++ ;
}
}
/* ---------------------------------------------------------------------- */
/* count the flops */
/* ---------------------------------------------------------------------- */
flop_count = 0.0 ;
for (k = 0 ; k < n ; k++)
{
/* off-diagonal nonzeros in column k of L: */
Lnz_k = Lp [k+1] - Lp [k] - 1 ;
Unz_k = Unz [k] - 1 ;
flop_count += (2 * Lnz_k * Unz_k) + Lnz_k ;
}
/* ---------------------------------------------------------------------- */
/* return the result */
/* ---------------------------------------------------------------------- */
plhs [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
pflop = mxGetPr (plhs [0]) ;
pflop [0] = flop_count ;
}
|
