-
Notifications
You must be signed in to change notification settings - Fork 1
/
compute_preconditioner.c
88 lines (67 loc) · 1.33 KB
/
compute_preconditioner.c
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
#include "header.h"
void compute_preconditioner(
mesh_struct *mesh,
matrix_int_crs_struct *splat_matrix,
int *blur_matrix[5],
double *Cn_matrix,
int *Cs_matrix,
double lambda,
double *con_arr,
double *M1_arr
)
/*
Compute M^-1, the (Jacobi) preconditioner
*/
{
arrayHeader *mesh_verts;
int maxInd;
int nbr_verts;
int vert;
double *diag_A;
double vert_con;
int ind;
int pixel;
double con;
/*
Get nbr of verts
*/
mesh_verts= mesh->mesh_verts;
maxInd= arrayMaxId(*mesh_verts);
nbr_verts= maxInd+1;
/*
Allocate memory for diag_A
*/
myCalloc(diag_A,double,nbr_verts,sizeof(double));
/*
Compute diag(A) =
lambda*(diag(Dm)-diag(Dn)*Bdiag*diag(Dn))+Sc
*/
for ( vert= 0 ; vert< nbr_verts ; vert++ ) {
diag_A[vert]= 0.0;
diag_A[vert]+= lambda*Cs_matrix[vert];
diag_A[vert]-= lambda*Cn_matrix[vert]*(2./4.)*Cn_matrix[vert];
vert_con= 0.0;
for ( ind= splat_matrix->row_ptr[vert] ;
ind< splat_matrix->row_ptr[vert+1] ;
ind++ ) {
pixel= splat_matrix->col_ind[ind];
con= con_arr[pixel];
vert_con+= con;
}
diag_A[vert]+= vert_con;
}
/*
Compute M^-1 =
diag(A)^-1
*/
for ( vert= 0 ; vert< nbr_verts ; vert++ ) {
if ( diag_A[vert] == 0 ) {
error_handler("pb !");
}
M1_arr[vert]= 1./diag_A[vert];
}
/*
Free diag_A
*/
myFree(diag_A);
}