[GRASS-dev] Re: [GRASS-user] v.lidar.edgedetection very slow
Jyothish Soman
jyothish.soman at gmail.com
Fri Jun 19 14:22:14 EDT 2009
Skipped content of type multipart/alternative-------------- next part --------------
#include <stdlib.h> /* imported libraries */
2 #include <stdio.h>
3 #include <math.h>
4 #include <grass/gis.h>
5 #include <grass/PolimiFunct.h>
6 #include<omp.h> /*Change here*/
7 /*--------------------------------------------------------------------------------------*/
8 /* Tcholetsky decomposition -> T= Lower Triangular Matrix */
9
10 void tcholDec(double **N, double **T, int n, int BW)
11 {
12 int i, j, k,num_p;
13 double somma;
14 omp_set_num_threads(omp_get_num_procs());
15 G_debug(3, "tcholDec(): n=%d BW=%d", n, BW);
16 #pragma omp for shared(i,N,BW,T) private(j,k,somma)
17 for (i = 0; i < n; i++) {
18 G_percent(i, n, 2);
19 for (j = 0; j < BW; j++) {
20 somma = N[i][j];
21 for (k = 1; k < BW; k++)
22 if (((i - k) >= 0) && ((j + k) < BW))
23 somma -= T[i - k][k] * T[i - k][j + k];
24 if (j == 0) {
25 if (somma <= 0.0)
26 G_fatal_error(_("Decomposition failed"));
27 T[i][0] = sqrt(somma);
28 }
29 else
30 T[i][j] = somma / T[i][0];
31 }
32 }
33
34 G_percent(i, n, 2);
35 return;
36 }
37
38 /*--------------------------------------------------------------------------------------*/
39 /* Tcholetsky matrix solution */
40
41 void tcholSolve(double **N, double *TN, double *parVect, int n, int BW)
42 {
43
44 double **T;
45 int i, j;
46 omp_set_num_threads(omp_get_num_procs());
47 /*--------------------------------------*/
48 T = G_alloc_matrix(n, BW);
49
50 /*--------------------------------------*/
51 tcholDec(N, T, n, BW); /* T computation */
52
53 /* Forward substitution */
54 parVect[0] = TN[0] / T[0][0];
#pragma omp for shared(i,TN,parvect) private(j)
55 for (i = 1; i < n; i++) {
56 parVect[i] = TN[i];
57 for (j = 0; j < i; j++)
58 if ((i - j) < BW)
59 parVect[i] -= T[j][i - j] * parVect[j];
60 parVect[i] = parVect[i] / T[i][0];
61 }
62 #pragma omp for shared(i) private()
63 /* Backward substitution */
64 parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
65 for (i = n - 2; i >= 0; i--) {
66 for (j = i + 1; j < n; j++)
67 if ((j - i) < BW)
68 parVect[i] -= T[i][j - i] * parVect[j];
69 parVect[i] = parVect[i] / T[i][0];
70 }
71
72 /*--------------------------------------*/
73 G_free_matrix(T);
74
75 return;
76 }
77
78
79 /*--------------------------------------------------------------------------------------*/
80 /* Soluzione con Tcholetsky -> la matrice T triangolare viene passata come paramtero e
81 non calcolata internamente alla procedura -> T = dmatrix (0, n-1, 0, BW-1) */
82
83 void tcholSolve2(double **N, double *TN, double **T, double *parVect, int n,
84 int BW)
85 {
86
87 int i, j;
88 omp_set_num_threads(omp_get_num_procs());
89 /* Forward substitution */
90 parVect[0] = TN[0] / T[0][0];
#pragma omp for shared(i) private()
91 for (i = 1; i < n; i++) {
92 parVect[i] = TN[i];
93 for (j = 0; j < i; j++)
94 if ((i - j) < BW)
95 parVect[i] -= T[j][i - j] * parVect[j];
96 parVect[i] = parVect[i] / T[i][0];
97 }
98
99 /* Backward substitution */
100 parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
#pragma omp for shared(i) private()
101 for (i = n - 2; i >= 0; i--) {
102 for (j = i + 1; j < n; j++)
103 if ((j - i) < BW)
104 parVect[i] -= T[i][j - i] * parVect[j];
105 parVect[i] = parVect[i] / T[i][0];
106 }
107
108 return;
109 }
110
111 /*--------------------------------------------------------------------------------------*/
112 /* Tcholetsky matrix invertion */
113
114 void tcholInv(double **N, double *invNdiag, int n, int BW)
115 {
116 double **T = NULL;
117 double *vect = NULL;
118 int i, j, k;
119 double somma;
120
121 /*--------------------------------------*/
122 G_alloc_matrix(n, BW);
123 G_alloc_vector(n);
124 omp_set_num_threads(omp_get_num_procs());
125 /* T computation */
126 tcholDec(N, T, n, BW);
127
128 /* T Diagonal invertion */
#pragma omp for shared(i)
129 for (i = 0; i < n; i++) {
130 T[i][0] = 1.0 / T[i][0];
131 }
132
133 /* N Diagonal invertion */
134 for (i = 0; i < n; i++) {
135 vect[0] = T[i][0];
136 invNdiag[i] = vect[0] * vect[0];
137 for (j = i + 1; j < n; j++) {
138 somma = 0.0;
139 for (k = i; k < j; k++) {
140 if ((j - k) < BW)
141 somma -= vect[k - i] * T[k][j - k];
142 }
143 vect[j - i] = somma * T[j][0];
144 invNdiag[i] += vect[j - i] * vect[j - i];
145 }
146 }
147
148 /*--------------------------------------*/
149 G_free_matrix(T);
150 G_free_vector(vect);
151
152 return;
153 }
154
155 /*--------------------------------------------------------------------------------------*/
156 /* Tcholetsky matrix solution and invertion */
157
158 void tcholSolveInv(double **N, double *TN, double *invNdiag, double *parVect,
159 int n, int BW)
160 {
161
162 double **T = NULL;
163 double *vect = NULL;
164 int i, j, k;
165 double somma;
166 omp_set_num_threads(omp_get_num_procs());
167 /*--------------------------------------*/
168 G_alloc_matrix(n, BW);
169 G_alloc_vector(n);
170
171 /* T computation */
172 tcholDec(N, T, n, BW);
173
174 /* Forward substitution */
175 parVect[0] = TN[0] / T[0][0];
#pragma omp for shared(i) private(j)
176 for (i = 1; i < n; i++) {
177 parVect[i] = TN[i];
178 for (j = 0; j < i; j++)
179 if ((i - j) < BW)
180 parVect[i] -= T[j][i - j] * parVect[j];
181 parVect[i] = parVect[i] / T[i][0];
182 }
183
184 /* Backward substitution */
185 parVect[n - 1] = parVect[n - 1] / T[n - 1][0];
#pragma omp for shared(i) private(j)
186 for (i = n - 2; i >= 0; i--) {
187 for (j = i + 1; j < n; j++)
188 if ((j - i) < BW)
189 parVect[i] -= T[i][j - i] * parVect[j];
190 parVect[i] = parVect[i] / T[i][0];
191 }
192
193 /* T Diagonal invertion */
#pragma omp for shared(i)
194 for (i = 0; i < n; i++) {
195 T[i][0] = 1.0 / T[i][0];
196 }
197
198 /* N Diagonal invertion */
199 for (i = 0; i < n; i++) {
200 vect[0] = T[i][0];
201 invNdiag[i] = vect[0] * vect[0];
202 for (j = i + 1; j < n; j++) {
203 somma = 0.0;
204 for (k = i; k < j; k++) {
205 if ((j - k) < BW)
206 somma -= vect[k - i] * T[k][j - k];
207 }
208 vect[j - i] = somma * T[j][0];
209 invNdiag[i] += vect[j - i] * vect[j - i];
210 }
211 }
212
213 /*--------------------------------------*/
214 G_free_matrix(T);
215 G_free_vector(vect);
216
217 return;
218 }
More information about the grass-dev
mailing list