IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_IKLU_Utils.h
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK_IKLU_UTILS_H
44 #define IFPACK_IKLU_UTILS_H
45 
46 #include <stdlib.h>
47 #include <limits.h>
48 #include <math.h>
49 #include <stdio.h>
50 
51 /* The code found in this file is adapted from CSparse Version 2.0.0
52  written by Tim Davis, UFL
53 */
54 
55 /* --- primary CSparse routines and data structures ------------------------- */
56 
57 typedef struct row_matrix /* matrix in compressed-row or triplet form */
58 {
59  int nzmax ; /* maximum number of entries */
60  int m ; /* number of rows */
61  int n ; /* number of columns */
62  int *p ; /* row pointers (size m+1) or col indices (size nzmax) */
63  int *j ; /* col indices, size nzmax */
64  double *x ; /* numerical values, size nzmax */
65  int nz ; /* # of entries in triplet matrix, -1 for compressed-row */
66 } csr ;
67 
68 csr *csr_add (const csr *A, const csr *B, double alpha, double beta) ;
69 csr *csr_multiply (const csr *A, const csr *B) ;
70 double csr_norm (const csr *A) ;
71 int csr_print (const csr *A, int brief) ;
72 csr *csr_transpose (const csr *A, int values) ;
73 
74 /* utilities */
75 void *csr_realloc (void *p, int n, size_t size, int *ok) ;
76 
77 /* csr utilities */
78 csr *csr_spalloc (int m, int n, int nzmax, int values, int triplet) ;
79 csr *csr_spfree (csr *A) ;
80 int csr_sprealloc (csr *A, int nzmax) ;
81 
82 
83 
84 /* --- secondary CSparse routines and data structures ----------------------- */
85 typedef struct cs_symbolic /* symbolic Cholesky, LU, or QR analysis */
86 {
87  int *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */
88  int *q ; /* fill-reducing column permutation for LU and QR */
89  int *parent ; /* elimination tree for Cholesky and QR */
90  int *cp ; /* column pointers for Cholesky, row counts for QR */
91  int *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */
92  int m2 ; /* # of rows for QR, after adding fictitious rows */
93  double lnz ; /* # entries in L for LU or Cholesky; in V for QR */
94  double unz ; /* # entries in U for LU; in R for QR */
95 } css ;
96 
97 typedef struct csr_numeric /* numeric Cholesky, LU, or QR factorization */
98 {
99  csr *L ; /* L for LU and Cholesky, V for QR */
100  csr *U ; /* U for LU, R for QR, not used for Cholesky */
101  int *pinv ; /* partial pivoting for LU */
102  int *perm ; /* partial pivoting for LU */
103  double *B ; /* beta [0..n-1] for QR */
104 } csrn ;
105 
106 typedef struct csr_dmperm_results /* csr_dmperm or csr_scc output */
107 {
108  int *p ; /* size m, row permutation */ /* may be back wards */
109  int *q ; /* size n, column permutation */
110  int *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */
111  int *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */
112  int nb ; /* # of blocks in fine dmperm decomposition */
113  int rr [5] ; /* coarse row decomposition */
114  int cc [5] ; /* coarse column decomposition */
115 } csrd ;
116 
117 int *csr_amd (int order, const csr *A) ;
118 
119 int csr_droptol (csr *A, double tol);
120 int csr_dropzeros (csr *A);
121 int csr_lsolve (const csr *L, double *x);
122 csrn *csr_lu (const csr *A, const css *S, double tol);
123 csr *csr_permute (const csr *A, const int *pinv, const int *q, int values);
124 css *csr_sqr (int order, const csr *A);
125 int csr_usolve (const csr *U, double *x);
126 
127 /* utilities */
128 css *csr_sfree (css *S) ;
129 csrd *csr_dfree (csrd *D);
130 csrn *csr_nfree (csrn *N);
131 
132 
133 
134 /* --- tertiary CSparse routines -------------------------------------------- */
135 double csr_cumsum (int *p, int *c, int n) ;
136 int csr_dfs (int j, csr *G, int top, int *xi, int *pstack, const int *pinv);
137 int csr_reach (csr *G, const csr *B, int k, int *xi, const int *pinv);
138 int csr_scatter (const csr *A, int j, double beta, int *w, double *x, int mark,
139  csr *C, int nz) ;
140 csrd *csr_scc (csr *A);
141 int csr_spsolve (csr *G, const csr *B, int k, int *xi,
142  double *x, const int *pinv, int up);
143 int csr_tdfs (int j, int k, int *head, const int *next, int *post,
144  int *stack) ;
145 /* utilities */
146 csrd *csr_dalloc (int m, int n);
147 /* csr utilities */
148 csrd *csr_ddone (csrd *D, csr *C, void *w, int ok) ;
149 csr *csr_done (csr *C, void *w, void *x, int ok) ;
150 int *csr_idone (int *p, csr *C, void *w, int ok) ;
151 csrn *csr_ndone (csrn *N, csr *C, void *w, void *x, int ok) ;
152 
153 int csr_fkeep (csr *A, int (*fkeep) (int, int, double, void *), void *other);
154 
155 #define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
156 #define CS_MIN(a,b) (((a) < (b)) ? (a) : (b))
157 #define CS_FLIP(i) (-(i)-2)
158 #define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
159 #define CS_MARKED(w,j) (w [j] < 0)
160 #define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; }
161 #define CS_CSC(A) (A && (A->nz == -1))
162 #define CS_TRIPLET(A) (A && (A->nz >= 0))
163 
164 #endif // IFPACK_IKLU_UTILS_H