Line data Source code
1 : #include "LG_internal.h"
2 : #include <float.h>
3 :
4 : // TODO: rename this method and add to src
5 : // TODO: add standard comments to the top of this file
6 : // TODO: what functions should be exposed to the user application?
7 :
8 : // candidates:
9 : // LAGraph_Laplacian: compute the Laplacian matrix
10 : // LAGraph_Happly: apply a Housefolder reflection
11 : // LAGraph_norm2: 2-norm of a vector
12 :
13 : // very specific to the HDIP method for computing the Fiedler vector:
14 : // LAGraph_hmhx: y = H*M*H*x = (M-u*x'-x*u)*x
15 : // LAGraph_mypcg2: preconditioned conjugate gradient (very specialized)
16 :
17 : //------------------------------------------------------------------------------
18 : // LAGraph_Hdip_Fiedler
19 : //------------------------------------------------------------------------------
20 :
21 : //-------------------------------------------------------------------------------------------------
22 : /*
23 : The happly function, applies a Householder Reflection
24 : y= happly(u,x,alpha)
25 : calculates y = H*x where H = (I - u*u'/alpha)
26 :
27 : H*x =(I - u*u'/alpha)*x=(x-u*(u'u*x)/alpha)
28 :
29 : y,u, and x are vectors of size n and alpha is a scalar.
30 : The inputs are all modified so y must be a different vector.
31 :
32 : Example:
33 : I=[4.0,5.0,6.0,7.0]
34 : v=gb.Vector.from_list(I)
35 : J=[1.0,2.0,3.0,4.0]
36 : x=gb.Vector.from_list(J)
37 : A = 2.0
38 : k=happly(v,x,A)
39 : print(k)
40 :
41 : This example, simply creates 2 vectors and uses a literal numerical value as alpha to calculate happly.
42 :
43 : Authors: Georgy Thayyil, Tim Davis, Mengting Cao
44 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
45 :
46 : */
47 722 : int LAGraph_Happly //happly Checked for pointer issues
48 : (
49 : //outputs:
50 : GrB_Vector y, // y output of Householder reflection on x.
51 : //inputs:
52 : GrB_Vector u, // u, the vector used for application of householder
53 : GrB_Vector x, // x, the vector on which householder reflection is applied
54 : float alpha, // the scalar alpha used for application of householder
55 : //error msg
56 : char *msg
57 : )
58 : {
59 : #if LAGRAPH_SUITESPARSE
60 722 : float reduced = 0.0;
61 : // y = u .* x
62 722 : GRB_TRY (GrB_eWiseAdd (y,NULL, NULL, GrB_TIMES_FP32, u, x, NULL));
63 : // reduced = sum (y)
64 722 : GRB_TRY (GrB_reduce(&reduced, NULL,GrB_PLUS_MONOID_FP32,y,NULL));
65 : // y = (-reduced/alpha)*u
66 722 : GRB_TRY (GrB_apply(y,NULL,NULL,GrB_TIMES_FP32,-reduced/alpha,u,NULL));
67 : // y = x + y
68 722 : GRB_TRY (GrB_eWiseAdd(y,NULL,NULL,GrB_PLUS_FP32,x,y,NULL));
69 :
70 722 : return (GrB_SUCCESS) ;
71 : #else
72 : return (GrB_NOT_IMPLEMENTED) ;
73 : #endif
74 : }
75 : //-------------------------------------------------------------------------------------------------
76 : /*
77 : hmhx = a function used to compute y = H*M*H*x = (M-u*x'-x*u)*x where x and y have the same size n
78 :
79 : Example:
80 : I=[4.0,5.0,6.0,7.0]
81 : v=gb.Vector.from_list(I)
82 : x=gb.Vector.from_list(I)
83 : A = 2.0
84 : I=[3,3,3]
85 : J=[1,2,3]
86 : K=[4,5,6]
87 : M=gb.Matrix.from_lists(I,J,K)
88 : M=M+M.transpose()
89 : mdiag=M.diag()
90 : indiag=1/mdiag.cast(gb.FP64)
91 : print(hmhx(indiag,v,x,A))
92 :
93 : This example starts off by creating 2 vectors and a matrix.
94 : Then makes sure the matrix is symmetric, and gets the inverse diagonal of the matrix.
95 : Lastly, the example prints, hmhx called using the modified matrix, and vectors.
96 :
97 :
98 : Authors: Georgy Thayyil, Tim Davis, Mengting Cao
99 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
100 :
101 : */
102 :
103 : #undef LG_FREE_WORK
104 : #define LG_FREE_WORK \
105 : { \
106 : /* free any workspace used here */ \
107 : GrB_free (&t) ; \
108 : }
109 :
110 : #undef LG_FREE_ALL
111 : #define LG_FREE_ALL \
112 : { \
113 : /* free any workspace used here */ \
114 : LG_FREE_WORK ; \
115 : /* free all the output variable(s) */ \
116 : /* take any other corrective action */ \
117 : }
118 :
119 361 : int LAGraph_hmhx //hmhx checked for pointer issues
120 : (
121 : //outputs:
122 : GrB_Vector z, //z output of hmhx
123 : //inputs:
124 : GrB_Matrix M, //Matrix used in hmhx
125 : GrB_Vector u, // Vector u used for happly
126 : GrB_Vector x, // Vector x used for happly
127 : float alpha, // the scalar alpha used for happly
128 : char *msg
129 : )
130 : {
131 : #if LAGRAPH_SUITESPARSE
132 361 : GrB_Vector t = NULL ;
133 361 : GrB_Index n = 0 ;
134 :
135 : // z = happly(u,x,alpha)
136 361 : LG_TRY (LAGraph_Happly(z,u,x,alpha,msg));
137 :
138 : // t = M*z
139 361 : GRB_TRY (GrB_Vector_size (&n, x)) ;
140 361 : GRB_TRY (GrB_Vector_new (&t, GrB_FP32, n)) ;
141 361 : GRB_TRY (GrB_mxv(t,NULL,NULL,GrB_PLUS_TIMES_SEMIRING_FP32,M,z,NULL)); // matrix vector multiply with z and M
142 :
143 : // z = happly (u,t,alpha)
144 361 : LG_TRY (LAGraph_Happly(z,u,t,alpha,msg));// z = happly with z,u and alpha
145 :
146 : // z(0) = 0
147 361 : GRB_TRY (GrB_Vector_setElement_FP32(z, 0, 0));
148 :
149 : // free workspace
150 361 : LG_FREE_WORK ;
151 361 : return (GrB_SUCCESS) ;
152 : #else
153 : return (GrB_NOT_IMPLEMENTED) ;
154 : #endif
155 : }
156 : //-------------------------------------------------------------------------------------------------
157 : /*
158 : norm2: This function aims to get the 2 norm of a vetor, and requires an import from python's math module.
159 :
160 : Example:
161 : I=[4.0,5.0,6.0,7.0]
162 : v=gb.Vector.from_list(I)
163 : k=norm2(v)
164 : print(k)
165 :
166 : Creates a vector, and passes it to the function and gets the norm2 of the function as a result.
167 :
168 : Authors: Georgy Thayyil, Tim Davis, Mengting Cao
169 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
170 : */
171 :
172 :
173 188 : int LAGraph_norm2 //norm2 checked for pointer mistakes
174 : (
175 : //outputs:
176 : float *norm2_helper,
177 : //inputs:
178 : GrB_Vector v,
179 : //error msg
180 : char *msg
181 :
182 : )
183 : {
184 : #if LAGRAPH_SUITESPARSE
185 188 : GrB_Vector t = NULL ;
186 : GrB_Index len ;
187 : float norm2;
188 :
189 188 : GRB_TRY (GrB_Vector_size (&len, v)) ;
190 188 : GRB_TRY (GrB_Vector_new (&t, GrB_FP32, len)) ;
191 :
192 : // t = v.^2
193 : // TODO: use GrB_TIMES_FP32 with GrB_eWiseMult for vanilla case
194 188 : GRB_TRY (GrB_apply (t, NULL, NULL, GxB_POW_FP32, v, (float) 2, NULL)) ;
195 :
196 188 : GRB_TRY (GrB_reduce (&norm2,NULL, GrB_PLUS_MONOID_FP32, t, NULL));
197 :
198 188 : norm2 = sqrtf (norm2) ;
199 :
200 188 : (*norm2_helper) = norm2 ;
201 188 : LG_FREE_WORK ;
202 188 : return (GrB_SUCCESS) ;
203 : #else
204 : return (GrB_NOT_IMPLEMENTED) ;
205 : #endif
206 : }
207 :
208 : //-------------------------------------------------------------------------------------------------
209 : /*
210 : Laplacian: This function computes the laplacian of a Matrix (This matrix has to be symmetric (a normal matrix, when added with its transpose becomes symmetric)), and returns a tuple with
211 : both the laplacian of the matrix, and the infinity norm of the matrix.
212 :
213 : Example:
214 : J=sorted(list(gb.Matrix.ssget(2760)), key=itemgetter(0))[0]
215 : h=J[1]+J[1].transpose()
216 : L=laplacian(h)
217 : print(L)
218 :
219 : The first line is simply the way to get a matrix from the SuiteSparse collection in python using the
220 : id of a matrix. It requires the import of "itemgetter" from operator, and the prior hand installation
221 : of ssgetpy.
222 :
223 : The second line simply makes sure that the matrix is symmetric by adding it with its transpose.
224 :
225 : The third line calls this laplacian function on the matrix and returns a tuple with the laplacian and the inf norm.
226 :
227 : Authors: Georgy Thayyil, Tim Davis, Mengting Cao
228 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
229 : */
230 :
231 : #undef LG_FREE_WORK
232 : #undef LG_FREE_ALL
233 : #define LG_FREE_WORK \
234 : { \
235 : /* free any workspace used here */ \
236 : GrB_free (&x); \
237 : GrB_free (&t); \
238 : GrB_free (&sparseM); \
239 : GrB_free (&DMatrix); \
240 : }
241 :
242 : #define LG_FREE_ALL \
243 : { \
244 : /* free any workspace used here */ \
245 : LG_FREE_WORK ; \
246 : /* free all the output variable(s) */ \
247 : GrB_free (&Lap) ; \
248 : }
249 :
250 3 : int LAGraph_Laplacian // compute the Laplacian matrix
251 : (
252 : // outputs:
253 : GrB_Matrix *Lap_handle, // the output Laplacian matrix
254 : float *infnorm, // infinity norm of Lap
255 : // inputs:
256 : GrB_Matrix G, // input matrix, symmetric
257 : // GrB_Type type, // the type of Lap, typically GrB_FP32, ...
258 : char *msg
259 : )
260 : {
261 : #if LAGRAPH_SUITESPARSE
262 : GrB_Index ncol;
263 3 : GrB_Vector x = NULL;
264 3 : GrB_Vector t = NULL;
265 3 : GrB_Matrix sparseM = NULL;
266 3 : GrB_Matrix DMatrix = NULL;
267 3 : GrB_Matrix Lap = NULL;
268 : // Assert G->A is symmetric
269 : // Lap = (float) offdiag (G->A)
270 3 : (*Lap_handle)= NULL ;
271 3 : GRB_TRY (GrB_Matrix_ncols(&ncol, G));
272 3 : GRB_TRY (GrB_Matrix_new (&Lap, GrB_FP32, ncol, ncol));
273 3 : GRB_TRY (GrB_select (Lap,NULL,NULL,GrB_OFFDIAG,G,0,NULL));
274 :
275 : // t = row degree of Lap
276 :
277 : // t = Lap * x via the LAGraph_plus_one_fp32 semiring
278 3 : GRB_TRY (GrB_Vector_new (&t, GrB_FP32, ncol));
279 3 : GRB_TRY (GrB_Vector_new (&x, GrB_FP32, ncol)) ;
280 : // create a dense vector filled with 0s
281 3 : GRB_TRY (GrB_assign (x, NULL, NULL, 0, GrB_ALL, ncol, NULL)) ;
282 : // x += Lap*x
283 3 : GRB_TRY(GrB_mxv(t,NULL,NULL,LAGraph_plus_one_fp32,Lap,x,NULL));
284 :
285 : //creates a sparse Matrix with same dimensions as Lap, and assigns -1 with Lap as a Mask
286 3 : GRB_TRY (GrB_Matrix_new (&sparseM, GrB_FP32, ncol, ncol));
287 : //Python code has descriptor = S indicating structural mask
288 3 : GRB_TRY(GrB_assign(sparseM, Lap, NULL,((double) -1),
289 : GrB_ALL, ncol, GrB_ALL, ncol, GrB_DESC_S));
290 :
291 : //create a mask of 0s in vector t, and use that to replace the 0s with 1s.
292 : //Python code has descriptor = S indicating structural mask
293 3 : GRB_TRY(GrB_assign(t,t,NULL,((double) 1),GrB_ALL,ncol,GrB_DESC_SC));
294 :
295 : //inf norm calc using vector d and MAX_MONOID
296 : float result ;
297 3 : GRB_TRY (GrB_reduce (&result, NULL, GrB_MAX_MONOID_FP32, t, NULL));
298 3 : result = result * 2 ;
299 :
300 : //Using Matrix_diag to create a diagonal matrix from a vector
301 3 : GRB_TRY (GrB_Matrix_diag(&DMatrix,t,0));
302 :
303 : //Calculating the Laplacian by adding the Dmatrix with SparseM.
304 3 : GRB_TRY (GrB_eWiseAdd (Lap, NULL, NULL, GrB_PLUS_FP32, DMatrix, sparseM, NULL));
305 :
306 : // free workspace and return result
307 3 : LG_FREE_WORK;
308 3 : (*Lap_handle)= Lap;
309 3 : (*infnorm)= result ;
310 3 : return (GrB_SUCCESS);
311 : #else
312 : return (GrB_NOT_IMPLEMENTED) ;
313 : #endif
314 : }
315 :
316 : //-------------------------------------------------------------------------------------------------
317 : /*
318 : mypcg2: Preconditioned conjugate gradient
319 :
320 : [x,k] = mypcg2 (L,u,v,alpha,C,b,tol,maxit)
321 :
322 : Solves A*x=b using pre-conditioned conjugage gradient and a diagonal
323 : preconditioner, where A = H*L*H and H is a Householder reflection with
324 : H = I - u*u'/alpha. A = H*L*H is equal to L-u*v'-v*u' but this is not
325 : formed explicitly. C is the pre-conditioner, and is a diagonal matrix with
326 : C(i,i) = 1/L(i,i). Thus, H*C*H is an approximation of the "inverse" A.
327 : More precisely if F = H*C*H then F is an approximation inv (A (2:n,2:n)),
328 : because A itself is exactly singular.
329 :
330 : L is the Laplacian matrix of an undirected graph, with L(i,i) = the degree of
331 : node i, and L(i,j) = -1 for each edge (i,j). L(i,i) must be > 0 for all i.
332 :
333 : Returns k as the # of iterations taken.
334 :
335 : From Golub and Van Loan, 2nd Edition, Algorithm 10.3.1, p529,
336 :
337 : Note that in the paper by Wu et al., the system A2*x=b has dimension n-1,
338 : with A2 = A (2:n,2:n). Here, all of A is handled, but A (:,1) and A (1,:)
339 : are all zero, and so is b (1) and the solution x (1). A2 must be symmetric
340 : and positive definite.
341 :
342 : In the GraphBLAS version, the input b can be overwritten with the solution x.
343 :
344 : Example:
345 : from operator import itemgetter
346 : J=sorted(list(gb.Matrix.ssget(2760)), key=itemgetter(0))[0]
347 : J[1].print(level=3)
348 : pattern=J[1].pattern(gb.types.FP32)
349 : h=pattern+pattern.transpose()
350 : omega=lap.laplacian(pattern)
351 : n=omega[0].nrows
352 : u=gb.Vector.dense(gb.types.FP32,n,fill=1)
353 : u[0]=1+sqrt(n)
354 : alpha=n+sqrt(n)
355 : c=omega[0].diag()
356 : indiag=(1/c.cast(gb.FP32)).diag()
357 : x=gb.Vector.dense(gb.types.FP32,n,fill=1)
358 : x[0]=0
359 : mypcg2=(omega[0],u,alpha,indiag,x,.000001,50)
360 : print(mypcg2)
361 :
362 : First we load in the matrix into python from suitesparse.
363 : Second we make sure that the matrix is of type FP32 (can be also be FP64)
364 : We make sure the matrix is symmetric, by adding it with its transpose
365 : Take the laplacian of that matrix.
366 : Create a vector filled with ones with the same size and number of rows in matrix
367 : Change first value in vector to (1+sqrt(num of rows in matrix))
368 : Create alpha, by (num of rows in matrix+sqrt(num of rows in matrix))
369 : Extract the diagonal of the modified matrix
370 : Create the inverse of that diagonal
371 : Create another vector filled with ones, with zero being its first value.
372 : Lastly call mypcg2 with the created variables as its parameters..
373 :
374 : Authors: Georgy Thayyil, Tim Davis, Mengting Cao
375 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
376 : */
377 :
378 : #undef LG_FREE_WORK
379 : #undef LG_FREE_ALL
380 : #define LG_FREE_WORK \
381 : { \
382 : /* free any workspace used here */ \
383 : GrB_free (&r); \
384 : GrB_free (&z); \
385 : GrB_free (&rho_helper); \
386 : GrB_free (&p); \
387 : GrB_free (&q); \
388 : GrB_free (&gamma_helper); \
389 : }
390 :
391 : #define LG_FREE_ALL \
392 : { \
393 : /* free any workspace used here */ \
394 : LG_FREE_WORK ; \
395 : /* free all the output variable(s) */ \
396 : GrB_free (&steper) ; \
397 : }
398 :
399 15 : int LAGraph_mypcg2
400 : (
401 : //outputs
402 : GrB_Vector *steper_handle, // TODO: do not create this; see usage below
403 : GrB_Index *k_result,
404 : // inputs:
405 : GrB_Matrix L, // input matrix, symmetric, result from Laplacian
406 : GrB_Vector u, //vector u will be passed into another function to create Householder reflection
407 : float malpha, //This float
408 : GrB_Matrix invdiag,
409 : GrB_Vector b,
410 : float tol,
411 : float maxit,
412 : //error msging
413 : char *msg
414 : )
415 : {
416 : #if LAGRAPH_SUITESPARSE
417 :
418 15 : (*steper_handle) = NULL ;
419 15 : GrB_Vector r = NULL ; // This vector will be a copy of vector b to make sure vector b remains unchanged.
420 15 : GrB_Vector z = NULL ; //used to apply preconditioner
421 15 : GrB_Vector rho_helper = NULL ; // used to help calculate rho
422 15 : GrB_Vector p = NULL ; //search direction
423 15 : GrB_Vector q = NULL ; //used for hmhx after finding next step in direction p
424 15 : GrB_Vector gamma_helper = NULL ;
425 15 : GrB_Vector steper = NULL ;
426 : GrB_Index n, k ;
427 : GrB_Index bsize;
428 : float rho;
429 : float rho_prior;
430 : float gamma;
431 : float stepsize;
432 : float rnorm;
433 : float beta;
434 :
435 : //Set n to be number of rows in Laplacian matrix
436 15 : GRB_TRY (GrB_Matrix_nrows(&n, L));
437 :
438 : //Set bsize to number of entries in vector b
439 15 : GRB_TRY (GrB_Vector_size(&bsize, b));
440 :
441 : //b[0] = 0 is required for the input
442 15 : GRB_TRY (GrB_Vector_setElement_FP32(b, 0, 0));
443 :
444 : //Set r to be equal to vector b
445 15 : GRB_TRY (GrB_Vector_dup(&r, b));
446 :
447 : //Set steper to size n, and filled with 0s.
448 15 : GRB_TRY (GrB_Vector_new (&steper, GrB_FP32, n));
449 : //GRB_TRY (GrB_apply (steper, NULL, NULL, GrB_FIRST_FP32, 0, steper,NULL));
450 15 : GRB_TRY(GrB_assign (steper, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
451 :
452 : //Initial rho
453 15 : rho = 1;
454 :
455 : //-------------''' Definition for helper vectors '''---------------------------
456 : //define rho_helper
457 15 : GRB_TRY (GrB_Vector_new (&rho_helper, GrB_FP32, bsize));
458 :
459 : //define the vector used to apply preconditioner
460 15 : GRB_TRY (GrB_Vector_new(&z, GrB_FP32, bsize));
461 :
462 : //define vector q used to do hmhx after next step in direction p
463 15 : GRB_TRY (GrB_Vector_new (&q, GrB_FP32, bsize));
464 :
465 : //define vector gamma_helper
466 15 : GRB_TRY (GrB_Vector_new (&gamma_helper, GrB_FP32, bsize));
467 :
468 : //------------------------------------------------------------------------------
469 :
470 173 : for (k=1;k<=maxit; k++)
471 : {
472 : //Apply the preconditioner, using hmhx
473 173 : LG_TRY (LAGraph_hmhx(z,invdiag,u,r,malpha,msg));
474 173 : GRB_TRY (GrB_Vector_setElement_FP32(z, 0, 0));
475 :
476 : //save the prior rho
477 173 : rho_prior=rho;
478 :
479 : //rho = r.emult(z).reduce_float(mon=gb.types.FP32.PLUS_MONOID) compute new rho
480 173 : GRB_TRY (GrB_eWiseAdd (rho_helper, NULL, NULL, GrB_TIMES_FP32, r, z, NULL));
481 173 : GRB_TRY (GrB_reduce(&rho,NULL,GrB_PLUS_MONOID_FP32,rho_helper,NULL));
482 :
483 173 : if(k==1){
484 : //first step is the direction p=z
485 15 : GRB_TRY (GrB_Vector_dup(&p, z));
486 : }else{
487 : //subsequent step in direction p
488 158 : beta = rho/rho_prior;
489 : //p=beta*p
490 158 : GRB_TRY (GrB_apply (p, NULL, NULL, GrB_TIMES_FP32, p, beta, NULL)) ;
491 : //p=p+z
492 158 : GRB_TRY (GrB_eWiseAdd(p, NULL, NULL, GrB_PLUS_FP32, p, z, NULL)) ;
493 : }
494 :
495 173 : GRB_TRY (GrB_Vector_setElement_FP32(p, 0, 0));
496 :
497 : // apply the matrix q = A*p
498 : //hmhx is used on q
499 173 : LG_TRY (LAGraph_hmhx(q,L,u,p,malpha,msg));
500 173 : GRB_TRY (GrB_Vector_setElement_FP32(q, 0, 0));
501 :
502 :
503 : //gamma=p.emult(q).reduce_float(mon=gb.types.FP32.PLUS_MONOID)
504 173 : GRB_TRY (GrB_eWiseMult (gamma_helper, NULL, NULL, GrB_TIMES_FP32, p, q, NULL));
505 173 : GRB_TRY (GrB_reduce(&gamma,NULL,GrB_PLUS_MONOID_FP32,gamma_helper,NULL));
506 :
507 173 : stepsize = rho/gamma;
508 :
509 : //takes a step towards the solution , steper = steper + stepsize*p
510 173 : GRB_TRY (GrB_apply (steper, NULL, GrB_PLUS_FP32, GrB_TIMES_FP32, p, stepsize, NULL));
511 :
512 : //r=r-stepsize*q
513 173 : GRB_TRY (GrB_apply (r, NULL, GrB_PLUS_FP32, GrB_TIMES_FP32, -stepsize, q, NULL));
514 :
515 173 : GRB_TRY (GrB_Vector_setElement_FP32(steper, 0, 0));
516 173 : GRB_TRY (GrB_Vector_setElement_FP32(r, 0, 0));
517 :
518 173 : LG_TRY (LAGraph_norm2(&rnorm,r,msg));// z = happly with z,u and alpha
519 :
520 173 : if(rnorm < tol){
521 15 : break;
522 : }
523 : }
524 :
525 : // free workspace and return result
526 15 : LG_FREE_WORK ;
527 15 : (*k_result) = k ;
528 15 : (*steper_handle) = steper ;
529 15 : return (GrB_SUCCESS);
530 : #else
531 : return (GrB_NOT_IMPLEMENTED) ;
532 : #endif
533 : }
534 :
535 :
536 : //-------------------------------------------------------------------------------------------------
537 : /*
538 : Hdip_fiedler: This function computes the Fiedler Vector of a graph by using the HDIP methos.
539 : x= hdip_fiedler(L) computes the Fiedler vector of the laplacian of a symmetric matrix.
540 :
541 : However, the laplacian function used here (created by the same authors) returns the
542 : laplacian of the matrix and its inf norm, in a tuple, as such the first input parameter
543 : to hdip_fiedler is a tuple with the laplacian of a symmetric matrix and its inf norm.
544 :
545 :
546 : The HDIP method used, was presented in this paper:
547 : ----------------------------------------------------------------------
548 : Jian-ping Wu, Jun-qiang Song, Wei-min Zhang, An efficient and
549 : accurate method to compute the Fiedler vector based on Householder
550 : deflation and inverse power iteration, Journal of Computational and
551 : Applied Mathematics, Volume 269, 2014, Pages 101-108, ISSN 0377-0427,
552 : https://doi.org/10.1016/j.cam.2014.03.018.
553 :
554 : Abstract: The Fiedler vector of a graph plays a vital role in many
555 : applications. But it is usually very expensive in that it involves
556 : the solution of an eigenvalue problem. In this paper, we introduce
557 : the inverse power method incorporated with the Householder deflation
558 : to compute the Fiedler Vector. In the inverse power iterations, the
559 : coefficient matrix is formed implicitly, to take advantage of the
560 : sparsity. The linear systems encountered at each iteration must be
561 : symmetric positive definite, thus the conjugate gradient method is
562 : used. In addition, preconditioning techniques are introduced to
563 : reduce the computation cost. Any kind of preconditioning techniques
564 : with dropping can be used. For the graphs related to some of the
565 : sparse matrices downloaded from the UF Sparse Matrix Collection, the
566 : experiments are compared to the known novel schemes. The results show
567 : that the provided method is more accurate. While it is slower than
568 : MC73 sequentially, it has good parallel efficiency compared with
569 : TraceMin. In addition, it is insensitive to the selection of
570 : parameters, which is superior to the other two methods.
571 : ----------------------------------------------------------------------
572 :
573 : Example:
574 : from operator import itemgetter
575 : J=sorted(list(gb.Matrix.ssget(2760)), key=itemgetter(0))[0]
576 : J[1].print(level=3)
577 : pattern=J[1].pattern(gb.types.FP32)
578 : h=pattern+pattern.transpose()
579 : omega=lap.laplacian(pattern)
580 : #omega.print(level=3)
581 : hdip2=hdip_fiedler(omega[0],omega[1])
582 : print(hdip2)
583 :
584 : The first line is simply the way to get a matrix from the SuiteSparse collection in python using the
585 : id of a matrix. It requires the import of "itemgetter" from operator, and the prior hand installation
586 :
587 : hdip2 from the example contains the vector as the first element in the tuple, the lambda value as the second element
588 : with the number of iterations as a tuple with a tuple in the third element.
589 :
590 : Authors: Georgy Thayyil, Tim Davis
591 : Acknowledgements: Michel Pelletier provided us with many helpful suggestions and assistance while developing this algorithm
592 : */
593 :
594 : #undef LG_FREE_WORK
595 : #undef LG_FREE_ALL
596 : #define LG_FREE_WORK \
597 : { \
598 : /* free any workspace used here */ \
599 : GrB_free (&u); \
600 : GrB_free (&y); \
601 : GrB_free (&lambhelper); \
602 : GrB_free (&indiag); \
603 : GrB_free (&x2) ; \
604 : }
605 :
606 : #define LG_FREE_ALL \
607 : { \
608 : /* free any workspace used here */ \
609 : LG_FREE_WORK ; \
610 : /* free all the output variable(s) */ \
611 : GrB_free (&x) ; \
612 : GrB_free (&iters) ; \
613 : }
614 :
615 : // TODO: need a basic API and an Advanced API
616 : // TODO: the basic API needs to use an LAGraph_Graph
617 : // TODO: the advance API can work on the Laplacian L and its infinity norm
618 : // TODO: this method has many inputs; the basic method should use defaults
619 : // TODO: need to add error checking of inputs
620 :
621 3 : int LAGraph_Hdip_Fiedler // compute the Hdip_Fiedler
622 : (
623 : // outputs:
624 : GrB_Vector *iters_handle, //This is a vector with number of inner and outer iterations listed in that order.
625 : float *lambda_result, // Lambda of hdip_fiedler
626 : GrB_Vector *x_handle, // the hdip fielder result vector
627 : // inputs:
628 : GrB_Matrix L, // input matrix, symmetric, result from Laplacian
629 : float InfNorm, // passed in from laplacian output
630 : GrB_Vector kmax, // a vector [20,50]
631 : float emax, // emax = .000001
632 : float tol, // tol = .000001
633 : // GrB_Type type, // the type of Lap, typically GrB_FP32, ...
634 : char *msg
635 : )
636 : {
637 : #if LAGRAPH_SUITESPARSE
638 :
639 : GrB_Index n;
640 3 : GrB_Index k_inner = 0 ;
641 : GrB_Index k_outer;
642 : GrB_Index kk; //used to hold output from mypcg2
643 : GrB_Index i; // This is the integer used in for loop
644 3 : GrB_Vector u = NULL;
645 3 : GrB_Vector y = NULL;
646 3 : GrB_Vector x = NULL, x2 = NULL ;
647 3 : GrB_Vector iters = NULL;
648 3 : GrB_Vector lambhelper = NULL;
649 : float alpha, lambda ;
650 : float last_err;
651 : float beta;
652 : float e;
653 : int kmaxZero; // kmax[0] set to 20
654 : int kmaxOne; // kmax[1] set to 50
655 3 : GrB_Matrix indiag = NULL ;
656 :
657 : //Set u(0) = 1+sqrt(n). u(1:n) = 1 and alpha = n+sqrt(n)
658 3 : GRB_TRY (GrB_Matrix_nrows(&n, L));
659 3 : GRB_TRY (GrB_Vector_new (&u, GrB_FP32, n));
660 : // u(0:n-1) = 1
661 3 : GRB_TRY (GrB_assign (u, NULL, NULL, 1, GrB_ALL, n, NULL)) ;
662 :
663 3 : GRB_TRY (GrB_Vector_setElement_FP32(u, 1+sqrtf(n), 0));
664 3 : GrB_wait (u, GrB_MATERIALIZE) ;
665 3 : alpha = n+sqrtf(n);
666 :
667 : //Set x(0) = 0 and x(1:n) = 1
668 3 : GRB_TRY (GrB_Vector_new (&x, GrB_FP32, n));
669 : // x(0:n-1) = 1
670 3 : GRB_TRY (GrB_assign (x, NULL, NULL, 1, GrB_ALL, n, NULL)) ;
671 3 : GRB_TRY (GrB_Vector_setElement_FP32(x, 0, 0));
672 :
673 : //indiag = diagonal matrix with indiag = 1/L[0], as a preconditioner
674 3 : GRB_TRY (GrB_Matrix_new(&indiag, GrB_FP32, n, n));
675 3 : GRB_TRY (GrB_select (indiag,NULL,NULL,GrB_DIAG, L, 0, NULL));
676 3 : GRB_TRY (GrB_apply (indiag, NULL, NULL, GrB_DIV_FP32, 1, indiag, NULL));
677 3 : last_err = FLT_MAX;
678 :
679 3 : GRB_TRY (GrB_Vector_new (&lambhelper, GrB_FP32, n));
680 :
681 : //used to set y = hmhx with m being L
682 3 : GRB_TRY (GrB_Vector_new (&y, GrB_FP32, n));
683 : //for i from 1 to kmax[0]+1
684 3 : GRB_TRY(GrB_Vector_extractElement(&kmaxZero,kmax,0));
685 : //setting up kmax[1]
686 3 : GRB_TRY(GrB_Vector_extractElement(&kmaxOne,kmax,1));
687 15 : for (i=1;i<=kmaxZero;i++)
688 : {
689 : //compute beta = 2-norm of x and set x = x/beta
690 15 : GRB_TRY (GrB_Vector_setElement_FP32(x, 0, 0));
691 15 : LG_TRY (LAGraph_norm2(&beta,x,msg));// z = happly with z,u and alpha
692 15 : GRB_TRY (GrB_apply (x, NULL, NULL, GrB_DIV_FP32, x,beta, NULL));
693 :
694 : //Set y = hmhx with m being L
695 15 : LG_TRY (LAGraph_hmhx(y,L,u,x,alpha,msg));
696 : //y(0)=0
697 15 : GRB_TRY (GrB_Vector_setElement_FP32(y, 0, 0));
698 :
699 : //lamb = x'*y
700 15 : GRB_TRY (GrB_eWiseMult (lambhelper, NULL, NULL, GrB_TIMES_FP32, x, y, NULL));
701 15 : GRB_TRY (GrB_reduce(&lambda,NULL,GrB_PLUS_MONOID_FP32,lambhelper,NULL));
702 :
703 : //getting the inf norm for the vector normer using norm(v,inf) = max(sum(abs(v))) vector v
704 : // e = norm(y-lambda*x, inf)
705 : //y = y - lambda*x
706 15 : GRB_TRY (GrB_apply (y, NULL, GrB_PLUS_FP32, GrB_TIMES_FP32, -lambda, x, NULL)) ;
707 : //getting abs(y)
708 15 : GRB_TRY (GrB_apply (lambhelper, NULL, NULL, GrB_ABS_FP32,y, NULL));
709 15 : GRB_TRY (GrB_reduce (&e, NULL, GrB_MAX_MONOID_FP32,y, NULL));
710 :
711 : //if e/inf norm of L<emax, break
712 15 : k_outer = i;
713 15 : e = e/InfNorm;
714 :
715 15 : if(e<emax || last_err<2*e)
716 : {
717 : break;
718 : }
719 12 : last_err=e;
720 :
721 : // TODO: revise x in place:
722 : //x=mypcg2(L,u,alpha,indiag,x,tol,kmax[1])
723 12 : LG_TRY (LAGraph_mypcg2(&x2,&kk,L,u,alpha,indiag,x,tol,kmaxOne,msg));
724 12 : GrB_free (&x) ;
725 12 : x = x2 ;
726 12 : x2 = NULL ;
727 12 : k_inner=k_inner+kk ;
728 :
729 12 : GRB_TRY (GrB_Vector_setElement_FP32(x, 0, 0));
730 :
731 : }
732 :
733 : //beta = u.emult(x).reduce_float(mon=gb.types.FP32.PLUS_MONOID)
734 3 : GRB_TRY (GrB_eWiseMult (lambhelper, NULL, NULL, GrB_TIMES_FP32, u,x, NULL));
735 3 : GRB_TRY (GrB_reduce(&beta,NULL,GrB_PLUS_MONOID_FP32,lambhelper,NULL));
736 3 : beta = beta/alpha;
737 :
738 : // set x=x-beta*u
739 3 : GRB_TRY (GrB_apply (x, NULL, GrB_PLUS_FP32, GrB_TIMES_FP32,-beta,u, NULL));
740 :
741 : //vectors start at 0
742 : //iters is used to return no. of inner and outer iterations
743 3 : GRB_TRY (GrB_Vector_new (&iters, GrB_FP32, 2));
744 3 : GRB_TRY (GrB_Vector_setElement_FP32(iters,k_outer,0));
745 3 : GRB_TRY (GrB_Vector_setElement_FP32(iters,k_inner,1));
746 :
747 : // free workspace and return result
748 3 : LG_FREE_WORK ;
749 3 : (*lambda_result) = lambda ;
750 3 : (*x_handle) = x;
751 3 : (*iters_handle) = iters ;
752 3 : return (GrB_SUCCESS);
753 : #else
754 : return (GrB_NOT_IMPLEMENTED) ;
755 : #endif
756 : }
757 :
|
