1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * Member functions for hypre_Vector class.
11 *
12 *****************************************************************************/
13
14 #include "_hypre_parcsr_mv.h"
15
16 HYPRE_Int hypre_FillResponseParToVectorAll(void*, HYPRE_Int, HYPRE_Int, void*, MPI_Comm, void**, HYPRE_Int*);
17
18 /*--------------------------------------------------------------------------
19 * hypre_ParVectorCreate
20 *--------------------------------------------------------------------------*/
21
22 /* If create is called and partitioning is NOT null, then it is assumed that it
23 is array of length 2 containing the start row of the calling processor
24 followed by the start row of the next processor - AHB 6/05 */
25
26 hypre_ParVector *
hypre_ParVectorCreate(MPI_Comm comm,HYPRE_BigInt global_size,HYPRE_BigInt * partitioning_in)27 hypre_ParVectorCreate( MPI_Comm comm,
28 HYPRE_BigInt global_size,
29 HYPRE_BigInt *partitioning_in )
30 {
31 hypre_ParVector *vector;
32 HYPRE_Int num_procs, my_id, local_size;
33 HYPRE_BigInt partitioning[2];
34
35 if (global_size < 0)
36 {
37 hypre_error_in_arg(2);
38 return NULL;
39 }
40 vector = hypre_CTAlloc(hypre_ParVector, 1, HYPRE_MEMORY_HOST);
41 hypre_MPI_Comm_rank(comm,&my_id);
42
43 if (!partitioning_in)
44 {
45 hypre_MPI_Comm_size(comm,&num_procs);
46 hypre_GenerateLocalPartitioning(global_size, num_procs, my_id, partitioning);
47 }
48 else
49 {
50 partitioning[0] = partitioning_in[0];
51 partitioning[1] = partitioning_in[1];
52 }
53 local_size = (HYPRE_Int) (partitioning[1] - partitioning[0]);
54
55 hypre_ParVectorAssumedPartition(vector) = NULL;
56
57 hypre_ParVectorComm(vector) = comm;
58 hypre_ParVectorGlobalSize(vector) = global_size;
59 hypre_ParVectorPartitioning(vector)[0] = partitioning[0];
60 hypre_ParVectorPartitioning(vector)[1] = partitioning[1];
61 hypre_ParVectorFirstIndex(vector) = hypre_ParVectorPartitioning(vector)[0];
62 hypre_ParVectorLastIndex(vector) = hypre_ParVectorPartitioning(vector)[1]-1;
63 hypre_ParVectorLocalVector(vector) = hypre_SeqVectorCreate(local_size);
64
65 /* set defaults */
66 hypre_ParVectorOwnsData(vector) = 1;
67 hypre_ParVectorActualLocalSize(vector) = 0;
68
69 return vector;
70 }
71
72 /*--------------------------------------------------------------------------
73 * hypre_ParMultiVectorCreate
74 *--------------------------------------------------------------------------*/
75
76 hypre_ParVector *
hypre_ParMultiVectorCreate(MPI_Comm comm,HYPRE_BigInt global_size,HYPRE_BigInt * partitioning,HYPRE_Int num_vectors)77 hypre_ParMultiVectorCreate( MPI_Comm comm,
78 HYPRE_BigInt global_size,
79 HYPRE_BigInt *partitioning,
80 HYPRE_Int num_vectors )
81 {
82 /* note that global_size is the global length of a single vector */
83 hypre_ParVector *vector = hypre_ParVectorCreate( comm, global_size, partitioning );
84 hypre_ParVectorNumVectors(vector) = num_vectors;
85 return vector;
86 }
87
88 /*--------------------------------------------------------------------------
89 * hypre_ParVectorDestroy
90 *--------------------------------------------------------------------------*/
91
92 HYPRE_Int
hypre_ParVectorDestroy(hypre_ParVector * vector)93 hypre_ParVectorDestroy( hypre_ParVector *vector )
94 {
95 if (vector)
96 {
97 if ( hypre_ParVectorOwnsData(vector) )
98 {
99 hypre_SeqVectorDestroy(hypre_ParVectorLocalVector(vector));
100 }
101
102 if (hypre_ParVectorAssumedPartition(vector))
103 {
104 hypre_AssumedPartitionDestroy(hypre_ParVectorAssumedPartition(vector));
105 }
106
107 hypre_TFree(vector, HYPRE_MEMORY_HOST);
108 }
109
110 return hypre_error_flag;
111 }
112
113 /*--------------------------------------------------------------------------
114 * hypre_ParVectorInitialize
115 *--------------------------------------------------------------------------*/
116
117 HYPRE_Int
hypre_ParVectorInitialize_v2(hypre_ParVector * vector,HYPRE_MemoryLocation memory_location)118 hypre_ParVectorInitialize_v2( hypre_ParVector *vector, HYPRE_MemoryLocation memory_location )
119 {
120 if (!vector)
121 {
122 hypre_error_in_arg(1);
123 return hypre_error_flag;
124 }
125 hypre_SeqVectorInitialize_v2(hypre_ParVectorLocalVector(vector), memory_location);
126
127 hypre_ParVectorActualLocalSize(vector) = hypre_VectorSize(hypre_ParVectorLocalVector(vector));
128
129 return hypre_error_flag;
130 }
131
132 HYPRE_Int
hypre_ParVectorInitialize(hypre_ParVector * vector)133 hypre_ParVectorInitialize( hypre_ParVector *vector )
134 {
135 return hypre_ParVectorInitialize_v2(vector, hypre_ParVectorMemoryLocation(vector));
136 }
137
138 /*--------------------------------------------------------------------------
139 * hypre_ParVectorSetDataOwner
140 *--------------------------------------------------------------------------*/
141
142 HYPRE_Int
hypre_ParVectorSetDataOwner(hypre_ParVector * vector,HYPRE_Int owns_data)143 hypre_ParVectorSetDataOwner( hypre_ParVector *vector,
144 HYPRE_Int owns_data )
145 {
146 if (!vector)
147 {
148 hypre_error_in_arg(1);
149 return hypre_error_flag;
150 }
151 hypre_ParVectorOwnsData(vector) = owns_data;
152
153 return hypre_error_flag;
154 }
155
156 /*--------------------------------------------------------------------------
157 * hypre_ParVectorSetNumVectors
158 * call before calling hypre_ParVectorInitialize
159 * probably this will do more harm than good, use hypre_ParMultiVectorCreate
160 *--------------------------------------------------------------------------*/
161 #if 0
162 HYPRE_Int
163 hypre_ParVectorSetNumVectors( hypre_ParVector *vector,
164 HYPRE_Int num_vectors )
165 {
166 HYPRE_Int ierr=0;
167 hypre_Vector *local_vector = hypre_ParVectorLocalVector(v);
168
169 hypre_SeqVectorSetNumVectors( local_vector, num_vectors );
170
171 return ierr;
172 }
173 #endif
174
175 /*--------------------------------------------------------------------------
176 * hypre_ParVectorRead
177 *--------------------------------------------------------------------------*/
178
179 hypre_ParVector*
hypre_ParVectorRead(MPI_Comm comm,const char * file_name)180 hypre_ParVectorRead( MPI_Comm comm,
181 const char *file_name )
182 {
183 char new_file_name[80];
184 hypre_ParVector *par_vector;
185 HYPRE_Int my_id;
186 HYPRE_BigInt partitioning[2];
187 HYPRE_BigInt global_size;
188 FILE *fp;
189
190 hypre_MPI_Comm_rank(comm, &my_id);
191
192 hypre_sprintf(new_file_name,"%s.INFO.%d",file_name,my_id);
193 fp = fopen(new_file_name, "r");
194 hypre_fscanf(fp, "%b\n", &global_size);
195 hypre_fscanf(fp, "%b\n", &partitioning[0]);
196 hypre_fscanf(fp, "%b\n", &partitioning[1]);
197 fclose (fp);
198 par_vector = hypre_CTAlloc(hypre_ParVector, 1, HYPRE_MEMORY_HOST);
199
200 hypre_ParVectorComm(par_vector) = comm;
201 hypre_ParVectorGlobalSize(par_vector) = global_size;
202
203 hypre_ParVectorFirstIndex(par_vector) = partitioning[0];
204 hypre_ParVectorLastIndex(par_vector) = partitioning[1]-1;
205
206 hypre_ParVectorPartitioning(par_vector)[0] = partitioning[0];
207 hypre_ParVectorPartitioning(par_vector)[1] = partitioning[1];
208
209 hypre_ParVectorOwnsData(par_vector) = 1;
210
211 hypre_sprintf(new_file_name,"%s.%d",file_name,my_id);
212 hypre_ParVectorLocalVector(par_vector) = hypre_SeqVectorRead(new_file_name);
213
214 /* multivector code not written yet */
215 hypre_assert( hypre_ParVectorNumVectors(par_vector) == 1 );
216
217 return par_vector;
218 }
219
220 /*--------------------------------------------------------------------------
221 * hypre_ParVectorPrint
222 *--------------------------------------------------------------------------*/
223
224 HYPRE_Int
hypre_ParVectorPrint(hypre_ParVector * vector,const char * file_name)225 hypre_ParVectorPrint( hypre_ParVector *vector,
226 const char *file_name )
227 {
228 char new_file_name[80];
229 hypre_Vector *local_vector;
230 MPI_Comm comm;
231 HYPRE_Int my_id;
232 HYPRE_BigInt *partitioning;
233 HYPRE_BigInt global_size;
234 FILE *fp;
235
236 if (!vector)
237 {
238 hypre_error_in_arg(1);
239 return hypre_error_flag;
240 }
241
242 local_vector = hypre_ParVectorLocalVector(vector);
243 comm = hypre_ParVectorComm(vector);
244 partitioning = hypre_ParVectorPartitioning(vector);
245 global_size = hypre_ParVectorGlobalSize(vector);
246
247 hypre_MPI_Comm_rank(comm,&my_id);
248 hypre_sprintf(new_file_name,"%s.%d",file_name,my_id);
249 hypre_SeqVectorPrint(local_vector,new_file_name);
250 hypre_sprintf(new_file_name,"%s.INFO.%d",file_name,my_id);
251 fp = fopen(new_file_name, "w");
252 hypre_fprintf(fp, "%b\n", global_size);
253 hypre_fprintf(fp, "%b\n", partitioning[0]);
254 hypre_fprintf(fp, "%b\n", partitioning[1]);
255
256 fclose(fp);
257
258 return hypre_error_flag;
259 }
260
261 /*--------------------------------------------------------------------------
262 * hypre_ParVectorSetConstantValues
263 *--------------------------------------------------------------------------*/
264
265 HYPRE_Int
hypre_ParVectorSetConstantValues(hypre_ParVector * v,HYPRE_Complex value)266 hypre_ParVectorSetConstantValues( hypre_ParVector *v,
267 HYPRE_Complex value )
268 {
269 hypre_Vector *v_local = hypre_ParVectorLocalVector(v);
270
271 return hypre_SeqVectorSetConstantValues(v_local,value);
272 }
273
274 /*--------------------------------------------------------------------------
275 * hypre_ParVectorSetRandomValues
276 *--------------------------------------------------------------------------*/
277
278 HYPRE_Int
hypre_ParVectorSetRandomValues(hypre_ParVector * v,HYPRE_Int seed)279 hypre_ParVectorSetRandomValues( hypre_ParVector *v,
280 HYPRE_Int seed )
281 {
282 HYPRE_Int my_id;
283 hypre_Vector *v_local = hypre_ParVectorLocalVector(v);
284
285 MPI_Comm comm = hypre_ParVectorComm(v);
286 hypre_MPI_Comm_rank(comm,&my_id);
287
288 seed *= (my_id+1);
289
290 return hypre_SeqVectorSetRandomValues(v_local, seed);
291 }
292
293 /*--------------------------------------------------------------------------
294 * hypre_ParVectorCopy
295 *--------------------------------------------------------------------------*/
296
297 HYPRE_Int
hypre_ParVectorCopy(hypre_ParVector * x,hypre_ParVector * y)298 hypre_ParVectorCopy( hypre_ParVector *x,
299 hypre_ParVector *y )
300 {
301 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
302 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
303 return hypre_SeqVectorCopy(x_local, y_local);
304 }
305
306 /*--------------------------------------------------------------------------
307 * hypre_ParVectorCloneShallow
308 * returns a complete copy of a hypre_ParVector x - a shallow copy, re-using
309 * the partitioning and data arrays of x
310 *--------------------------------------------------------------------------*/
311
312 hypre_ParVector *
hypre_ParVectorCloneShallow(hypre_ParVector * x)313 hypre_ParVectorCloneShallow( hypre_ParVector *x )
314 {
315 hypre_ParVector * y =
316 hypre_ParVectorCreate(hypre_ParVectorComm(x), hypre_ParVectorGlobalSize(x),
317 hypre_ParVectorPartitioning(x));
318
319 hypre_ParVectorOwnsData(y) = 1;
320 /* ...This vector owns its local vector, although the local vector doesn't
321 * own _its_ data */
322 hypre_SeqVectorDestroy( hypre_ParVectorLocalVector(y) );
323 hypre_ParVectorLocalVector(y) = hypre_SeqVectorCloneShallow(hypre_ParVectorLocalVector(x) );
324 hypre_ParVectorFirstIndex(y) = hypre_ParVectorFirstIndex(x);
325
326 return y;
327 }
328
329 hypre_ParVector *
hypre_ParVectorCloneDeep_v2(hypre_ParVector * x,HYPRE_MemoryLocation memory_location)330 hypre_ParVectorCloneDeep_v2( hypre_ParVector *x, HYPRE_MemoryLocation memory_location )
331 {
332 hypre_ParVector *y =
333 hypre_ParVectorCreate(hypre_ParVectorComm(x), hypre_ParVectorGlobalSize(x),
334 hypre_ParVectorPartitioning(x));
335
336 hypre_ParVectorOwnsData(y) = 1;
337 hypre_SeqVectorDestroy( hypre_ParVectorLocalVector(y) );
338 hypre_ParVectorLocalVector(y) = hypre_SeqVectorCloneDeep_v2( hypre_ParVectorLocalVector(x),
339 memory_location );
340 hypre_ParVectorFirstIndex(y) = hypre_ParVectorFirstIndex(x); //RL: WHY HERE?
341
342 return y;
343 }
344
345 HYPRE_Int
hypre_ParVectorMigrate(hypre_ParVector * x,HYPRE_MemoryLocation memory_location)346 hypre_ParVectorMigrate(hypre_ParVector *x, HYPRE_MemoryLocation memory_location)
347 {
348 if (!x)
349 {
350 return hypre_error_flag;
351 }
352
353 if ( hypre_GetActualMemLocation(memory_location) !=
354 hypre_GetActualMemLocation(hypre_ParVectorMemoryLocation(x)) )
355 {
356 hypre_Vector *x_local = hypre_SeqVectorCloneDeep_v2(hypre_ParVectorLocalVector(x), memory_location);
357 hypre_SeqVectorDestroy(hypre_ParVectorLocalVector(x));
358 hypre_ParVectorLocalVector(x) = x_local;
359 }
360 else
361 {
362 hypre_VectorMemoryLocation(hypre_ParVectorLocalVector(x)) = memory_location;
363 }
364
365 return hypre_error_flag;
366 }
367
368
369 /*--------------------------------------------------------------------------
370 * hypre_ParVectorScale
371 *--------------------------------------------------------------------------*/
372
373 HYPRE_Int
hypre_ParVectorScale(HYPRE_Complex alpha,hypre_ParVector * y)374 hypre_ParVectorScale( HYPRE_Complex alpha,
375 hypre_ParVector *y )
376 {
377 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
378
379 return hypre_SeqVectorScale( alpha, y_local);
380 }
381
382 /*--------------------------------------------------------------------------
383 * hypre_ParVectorAxpy
384 *--------------------------------------------------------------------------*/
385
386 HYPRE_Int
hypre_ParVectorAxpy(HYPRE_Complex alpha,hypre_ParVector * x,hypre_ParVector * y)387 hypre_ParVectorAxpy( HYPRE_Complex alpha,
388 hypre_ParVector *x,
389 hypre_ParVector *y )
390 {
391 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
392 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
393
394 return hypre_SeqVectorAxpy( alpha, x_local, y_local);
395 }
396
397 /*--------------------------------------------------------------------------
398 * hypre_ParVectorInnerProd
399 *--------------------------------------------------------------------------*/
400
401 HYPRE_Real
hypre_ParVectorInnerProd(hypre_ParVector * x,hypre_ParVector * y)402 hypre_ParVectorInnerProd( hypre_ParVector *x,
403 hypre_ParVector *y )
404 {
405 MPI_Comm comm = hypre_ParVectorComm(x);
406 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
407 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
408
409 HYPRE_Real result = 0.0;
410 HYPRE_Real local_result = hypre_SeqVectorInnerProd(x_local, y_local);
411
412 #ifdef HYPRE_PROFILE
413 hypre_profile_times[HYPRE_TIMER_ID_ALL_REDUCE] -= hypre_MPI_Wtime();
414 #endif
415 hypre_MPI_Allreduce(&local_result, &result, 1, HYPRE_MPI_REAL,
416 hypre_MPI_SUM, comm);
417 #ifdef HYPRE_PROFILE
418 hypre_profile_times[HYPRE_TIMER_ID_ALL_REDUCE] += hypre_MPI_Wtime();
419 #endif
420
421 return result;
422 }
423
424 /*--------------------------------------------------------------------------
425 * hypre_ParVectorElmdivpy
426 * y = y + x ./ b [MATLAB Notation]
427 *--------------------------------------------------------------------------*/
428
429 HYPRE_Int
hypre_ParVectorElmdivpy(hypre_ParVector * x,hypre_ParVector * b,hypre_ParVector * y)430 hypre_ParVectorElmdivpy( hypre_ParVector *x,
431 hypre_ParVector *b,
432 hypre_ParVector *y )
433 {
434 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
435 hypre_Vector *b_local = hypre_ParVectorLocalVector(b);
436 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
437
438 return hypre_SeqVectorElmdivpy(x_local, b_local, y_local);
439 }
440
441 /*--------------------------------------------------------------------------
442 * hypre_ParVectorElmdivpyMarked
443 * y[i] += x[i] / b[i] where marker[i] == marker_val
444 *--------------------------------------------------------------------------*/
445
446 HYPRE_Int
hypre_ParVectorElmdivpyMarked(hypre_ParVector * x,hypre_ParVector * b,hypre_ParVector * y,HYPRE_Int * marker,HYPRE_Int marker_val)447 hypre_ParVectorElmdivpyMarked( hypre_ParVector *x,
448 hypre_ParVector *b,
449 hypre_ParVector *y,
450 HYPRE_Int *marker,
451 HYPRE_Int marker_val )
452 {
453 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
454 hypre_Vector *b_local = hypre_ParVectorLocalVector(b);
455 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
456
457 return hypre_SeqVectorElmdivpyMarked(x_local, b_local, y_local, marker, marker_val);
458 }
459
460 /*--------------------------------------------------------------------------
461 * hypre_VectorToParVector:
462 * generates a ParVector from a Vector on proc 0 and distributes the pieces
463 * to the other procs in comm
464 *--------------------------------------------------------------------------*/
465
466 hypre_ParVector *
hypre_VectorToParVector(MPI_Comm comm,hypre_Vector * v,HYPRE_BigInt * vec_starts)467 hypre_VectorToParVector ( MPI_Comm comm,
468 hypre_Vector *v,
469 HYPRE_BigInt *vec_starts )
470 {
471 HYPRE_BigInt global_size;
472 HYPRE_BigInt *global_vec_starts = NULL;
473 HYPRE_BigInt first_index;
474 HYPRE_BigInt last_index;
475 HYPRE_Int local_size;
476 HYPRE_Int num_vectors;
477 HYPRE_Int num_procs, my_id;
478 HYPRE_Int global_vecstride, vecstride, idxstride;
479 hypre_ParVector *par_vector;
480 hypre_Vector *local_vector;
481 HYPRE_Complex *v_data;
482 HYPRE_Complex *local_data;
483 hypre_MPI_Request *requests;
484 hypre_MPI_Status *status, status0;
485 HYPRE_Int i, j, k, p;
486
487 hypre_MPI_Comm_size(comm, &num_procs);
488 hypre_MPI_Comm_rank(comm, &my_id);
489
490 if (my_id == 0)
491 {
492 global_size = (HYPRE_BigInt)hypre_VectorSize(v);
493 v_data = hypre_VectorData(v);
494 num_vectors = hypre_VectorNumVectors(v); /* for multivectors */
495 global_vecstride = hypre_VectorVectorStride(v);
496 }
497
498 hypre_MPI_Bcast(&global_size,1,HYPRE_MPI_INT,0,comm);
499 hypre_MPI_Bcast(&num_vectors,1,HYPRE_MPI_INT,0,comm);
500 hypre_MPI_Bcast(&global_vecstride,1,HYPRE_MPI_INT,0,comm);
501
502 if (num_vectors == 1)
503 {
504 par_vector = hypre_ParVectorCreate(comm, global_size, vec_starts);
505 }
506 else
507 {
508 par_vector = hypre_ParMultiVectorCreate(comm, global_size, vec_starts, num_vectors);
509 }
510
511 vec_starts = hypre_ParVectorPartitioning(par_vector);
512 first_index = hypre_ParVectorFirstIndex(par_vector);
513 last_index = hypre_ParVectorLastIndex(par_vector);
514 local_size = (HYPRE_Int)(last_index - first_index) + 1;
515
516 if (my_id == 0)
517 {
518 global_vec_starts = hypre_CTAlloc(HYPRE_BigInt, num_procs+1, HYPRE_MEMORY_HOST);
519 }
520 hypre_MPI_Gather(&first_index, 1, HYPRE_MPI_BIG_INT, global_vec_starts,
521 1, HYPRE_MPI_BIG_INT, 0, comm);
522 if (my_id == 0)
523 {
524 global_vec_starts[num_procs] = hypre_ParVectorGlobalSize(par_vector);
525 }
526
527 hypre_ParVectorInitialize(par_vector);
528 local_vector = hypre_ParVectorLocalVector(par_vector);
529 local_data = hypre_VectorData(local_vector);
530 vecstride = hypre_VectorVectorStride(local_vector);
531 idxstride = hypre_VectorIndexStride(local_vector);
532 /* so far the only implemented multivector StorageMethod is 0 */
533 hypre_assert( idxstride==1 );
534
535 if (my_id == 0)
536 {
537 requests = hypre_CTAlloc(hypre_MPI_Request, num_vectors*(num_procs-1), HYPRE_MEMORY_HOST);
538 status = hypre_CTAlloc(hypre_MPI_Status, num_vectors*(num_procs-1), HYPRE_MEMORY_HOST);
539 k = 0;
540 for (p = 1; p<num_procs; p++)
541 for (j = 0; j<num_vectors; ++j)
542 {
543 hypre_MPI_Isend( &v_data[(HYPRE_Int) global_vec_starts[p]] + j*global_vecstride,
544 (HYPRE_Int)(global_vec_starts[p+1] - global_vec_starts[p]),
545 HYPRE_MPI_COMPLEX, p, 0, comm, &requests[k++] );
546 }
547 if (num_vectors == 1)
548 {
549 for (i = 0; i < local_size; i++)
550 local_data[i] = v_data[i];
551 }
552 else
553 {
554 for (j = 0; j<num_vectors; ++j)
555 {
556 for (i = 0; i < local_size; i++)
557 local_data[i+j*vecstride] = v_data[i+j*global_vecstride];
558 }
559 }
560 hypre_MPI_Waitall(num_procs-1,requests, status);
561 hypre_TFree(requests, HYPRE_MEMORY_HOST);
562 hypre_TFree(status, HYPRE_MEMORY_HOST);
563 }
564 else
565 {
566 for ( j=0; j<num_vectors; ++j )
567 hypre_MPI_Recv( local_data+j*vecstride, local_size, HYPRE_MPI_COMPLEX,
568 0, 0, comm,&status0 );
569 }
570
571 if (global_vec_starts)
572 {
573 hypre_TFree(global_vec_starts, HYPRE_MEMORY_HOST);
574 }
575
576 return par_vector;
577 }
578
579 /*--------------------------------------------------------------------------
580 * hypre_ParVectorToVectorAll:
581 * generates a Vector on every proc which has a piece of the data
582 * from a ParVector on several procs in comm,
583 * vec_starts needs to contain the partitioning across all procs in comm
584 *--------------------------------------------------------------------------*/
585
586 hypre_Vector *
hypre_ParVectorToVectorAll(hypre_ParVector * par_v)587 hypre_ParVectorToVectorAll( hypre_ParVector *par_v )
588 {
589 MPI_Comm comm = hypre_ParVectorComm(par_v);
590 HYPRE_BigInt global_size = hypre_ParVectorGlobalSize(par_v);
591 hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_v);
592 HYPRE_Int num_procs, my_id;
593 HYPRE_Int num_vectors = hypre_ParVectorNumVectors(par_v);
594 hypre_Vector *vector;
595 HYPRE_Complex *vector_data;
596 HYPRE_Complex *local_data;
597 HYPRE_Int local_size;
598 hypre_MPI_Request *requests;
599 hypre_MPI_Status *status;
600 HYPRE_Int i, j;
601 HYPRE_Int *used_procs;
602 HYPRE_Int num_types, num_requests;
603 HYPRE_Int vec_len, proc_id;
604
605 HYPRE_Int *new_vec_starts;
606
607 HYPRE_Int num_contacts;
608 HYPRE_Int contact_proc_list[1];
609 HYPRE_Int contact_send_buf[1];
610 HYPRE_Int contact_send_buf_starts[2];
611 HYPRE_Int max_response_size;
612 HYPRE_Int *response_recv_buf=NULL;
613 HYPRE_Int *response_recv_buf_starts = NULL;
614 hypre_DataExchangeResponse response_obj;
615 hypre_ProcListElements send_proc_obj;
616
617 HYPRE_Int *send_info = NULL;
618 hypre_MPI_Status status1;
619 HYPRE_Int count, tag1 = 112, tag2 = 223;
620 HYPRE_Int start;
621
622 hypre_MPI_Comm_size(comm, &num_procs);
623 hypre_MPI_Comm_rank(comm, &my_id);
624
625 local_size = (HYPRE_Int)(hypre_ParVectorLastIndex(par_v) -
626 hypre_ParVectorFirstIndex(par_v) + 1);
627
628 /* determine procs which hold data of par_v and store ids in used_procs */
629 /* we need to do an exchange data for this. If I own row then I will contact
630 processor 0 with the endpoint of my local range */
631
632 if (local_size > 0)
633 {
634 num_contacts = 1;
635 contact_proc_list[0] = 0;
636 contact_send_buf[0] = hypre_ParVectorLastIndex(par_v);
637 contact_send_buf_starts[0] = 0;
638 contact_send_buf_starts[1] = 1;
639 }
640 else
641 {
642 num_contacts = 0;
643 contact_send_buf_starts[0] = 0;
644 contact_send_buf_starts[1] = 0;
645 }
646
647 /*build the response object*/
648 /*send_proc_obj will be for saving info from contacts */
649 send_proc_obj.length = 0;
650 send_proc_obj.storage_length = 10;
651 send_proc_obj.id = hypre_CTAlloc(HYPRE_Int, send_proc_obj.storage_length, HYPRE_MEMORY_HOST);
652 send_proc_obj.vec_starts =
653 hypre_CTAlloc(HYPRE_Int, send_proc_obj.storage_length + 1, HYPRE_MEMORY_HOST);
654 send_proc_obj.vec_starts[0] = 0;
655 send_proc_obj.element_storage_length = 10;
656 send_proc_obj.elements =
657 hypre_CTAlloc(HYPRE_BigInt, send_proc_obj.element_storage_length, HYPRE_MEMORY_HOST);
658
659 max_response_size = 0; /* each response is null */
660 response_obj.fill_response = hypre_FillResponseParToVectorAll;
661 response_obj.data1 = NULL;
662 response_obj.data2 = &send_proc_obj; /*this is where we keep info from contacts*/
663
664
665 hypre_DataExchangeList(num_contacts,
666 contact_proc_list, contact_send_buf,
667 contact_send_buf_starts, sizeof(HYPRE_Int),
668 //0, &response_obj,
669 sizeof(HYPRE_Int), &response_obj,
670 max_response_size, 1,
671 comm, (void**) &response_recv_buf,
672 &response_recv_buf_starts);
673
674 /* now processor 0 should have a list of ranges for processors that have rows -
675 these are in send_proc_obj - it needs to create the new list of processors
676 and also an array of vec starts - and send to those who own row*/
677 if (my_id)
678 {
679 if (local_size)
680 {
681 /* look for a message from processor 0 */
682 hypre_MPI_Probe(0, tag1, comm, &status1);
683 hypre_MPI_Get_count(&status1, HYPRE_MPI_INT, &count);
684
685 send_info = hypre_CTAlloc(HYPRE_Int, count, HYPRE_MEMORY_HOST);
686 hypre_MPI_Recv(send_info, count, HYPRE_MPI_INT, 0, tag1, comm, &status1);
687
688 /* now unpack */
689 num_types = send_info[0];
690 used_procs = hypre_CTAlloc(HYPRE_Int, num_types, HYPRE_MEMORY_HOST);
691 new_vec_starts = hypre_CTAlloc(HYPRE_Int, num_types+1, HYPRE_MEMORY_HOST);
692
693 for (i=1; i<= num_types; i++)
694 {
695 used_procs[i-1] = (HYPRE_Int)send_info[i];
696 }
697 for (i=num_types+1; i< count; i++)
698 {
699 new_vec_starts[i-num_types-1] = send_info[i] ;
700 }
701 }
702 else /* clean up and exit */
703 {
704 hypre_TFree(send_proc_obj.vec_starts, HYPRE_MEMORY_HOST);
705 hypre_TFree(send_proc_obj.id, HYPRE_MEMORY_HOST);
706 hypre_TFree(send_proc_obj.elements, HYPRE_MEMORY_HOST);
707 if(response_recv_buf) hypre_TFree(response_recv_buf, HYPRE_MEMORY_HOST);
708 if(response_recv_buf_starts) hypre_TFree(response_recv_buf_starts, HYPRE_MEMORY_HOST);
709 return NULL;
710 }
711 }
712 else /* my_id ==0 */
713 {
714 num_types = send_proc_obj.length;
715 used_procs = hypre_CTAlloc(HYPRE_Int, num_types, HYPRE_MEMORY_HOST);
716 new_vec_starts = hypre_CTAlloc(HYPRE_Int, num_types+1, HYPRE_MEMORY_HOST);
717
718 new_vec_starts[0] = 0;
719 for (i=0; i< num_types; i++)
720 {
721 used_procs[i] = send_proc_obj.id[i];
722 new_vec_starts[i+1] = send_proc_obj.elements[i]+1;
723 }
724 hypre_qsort0(used_procs, 0, num_types-1);
725 hypre_qsort0(new_vec_starts, 0, num_types);
726 /*now we need to put into an array to send */
727 count = 2*num_types+2;
728 send_info = hypre_CTAlloc(HYPRE_Int, count, HYPRE_MEMORY_HOST);
729 send_info[0] = num_types;
730 for (i=1; i<= num_types; i++)
731 {
732 send_info[i] = (HYPRE_Int)used_procs[i-1];
733 }
734 for (i=num_types+1; i< count; i++)
735 {
736 send_info[i] = new_vec_starts[i-num_types-1];
737 }
738 requests = hypre_CTAlloc(hypre_MPI_Request, num_types, HYPRE_MEMORY_HOST);
739 status = hypre_CTAlloc(hypre_MPI_Status, num_types, HYPRE_MEMORY_HOST);
740
741 /* don't send to myself - these are sorted so my id would be first*/
742 start = 0;
743 if (used_procs[0] == 0)
744 {
745 start = 1;
746 }
747
748
749 for (i=start; i < num_types; i++)
750 {
751 hypre_MPI_Isend(send_info, count, HYPRE_MPI_INT, used_procs[i],
752 tag1, comm, &requests[i-start]);
753 }
754 hypre_MPI_Waitall(num_types-start, requests, status);
755
756 hypre_TFree(status, HYPRE_MEMORY_HOST);
757 hypre_TFree(requests, HYPRE_MEMORY_HOST);
758 }
759
760 /* clean up */
761 hypre_TFree(send_proc_obj.vec_starts, HYPRE_MEMORY_HOST);
762 hypre_TFree(send_proc_obj.id, HYPRE_MEMORY_HOST);
763 hypre_TFree(send_proc_obj.elements, HYPRE_MEMORY_HOST);
764 hypre_TFree(send_info, HYPRE_MEMORY_HOST);
765 if(response_recv_buf) hypre_TFree(response_recv_buf, HYPRE_MEMORY_HOST);
766 if(response_recv_buf_starts) hypre_TFree(response_recv_buf_starts, HYPRE_MEMORY_HOST);
767
768 /* now proc 0 can exit if it has no rows */
769 if (!local_size) {
770 hypre_TFree(used_procs, HYPRE_MEMORY_HOST);
771 hypre_TFree(new_vec_starts, HYPRE_MEMORY_HOST);
772 return NULL;
773 }
774
775 /* everyone left has rows and knows: new_vec_starts, num_types, and used_procs */
776
777 /* this vector should be rather small */
778
779 local_data = hypre_VectorData(local_vector);
780 vector = hypre_SeqVectorCreate((HYPRE_Int)global_size);
781 hypre_VectorNumVectors(vector) = num_vectors;
782 hypre_SeqVectorInitialize(vector);
783 vector_data = hypre_VectorData(vector);
784
785 num_requests = 2*num_types;
786
787 requests = hypre_CTAlloc(hypre_MPI_Request, num_requests, HYPRE_MEMORY_HOST);
788 status = hypre_CTAlloc(hypre_MPI_Status, num_requests, HYPRE_MEMORY_HOST);
789
790 /* initialize data exchange among used_procs and generate vector - here we
791 send to ourself also*/
792
793 j = 0;
794 for (i = 0; i < num_types; i++)
795 {
796 proc_id = used_procs[i];
797 vec_len = (HYPRE_Int)(new_vec_starts[i+1] - new_vec_starts[i]);
798 hypre_MPI_Irecv(&vector_data[(HYPRE_Int)new_vec_starts[i]], num_vectors*vec_len,
799 HYPRE_MPI_COMPLEX, proc_id, tag2, comm, &requests[j++]);
800 }
801 for (i = 0; i < num_types; i++)
802 {
803 hypre_MPI_Isend(local_data, num_vectors*local_size, HYPRE_MPI_COMPLEX,
804 used_procs[i], tag2, comm, &requests[j++]);
805 }
806
807 hypre_MPI_Waitall(num_requests, requests, status);
808
809 if (num_requests)
810 {
811 hypre_TFree(requests, HYPRE_MEMORY_HOST);
812 hypre_TFree(status, HYPRE_MEMORY_HOST);
813 hypre_TFree(used_procs, HYPRE_MEMORY_HOST);
814 }
815
816 hypre_TFree(new_vec_starts, HYPRE_MEMORY_HOST);
817
818 return vector;
819 }
820
821 /*--------------------------------------------------------------------------
822 * hypre_ParVectorPrintIJ
823 *--------------------------------------------------------------------------*/
824
825 HYPRE_Int
hypre_ParVectorPrintIJ(hypre_ParVector * vector,HYPRE_Int base_j,const char * filename)826 hypre_ParVectorPrintIJ( hypre_ParVector *vector,
827 HYPRE_Int base_j,
828 const char *filename )
829 {
830 MPI_Comm comm;
831 HYPRE_BigInt global_size, j;
832 HYPRE_BigInt *partitioning;
833 HYPRE_Complex *local_data;
834 HYPRE_Int myid, num_procs, i, part0;
835 char new_filename[255];
836 FILE *file;
837 if (!vector)
838 {
839 hypre_error_in_arg(1);
840 return hypre_error_flag;
841 }
842 comm = hypre_ParVectorComm(vector);
843 global_size = hypre_ParVectorGlobalSize(vector);
844 partitioning = hypre_ParVectorPartitioning(vector);
845
846 /* multivector code not written yet */
847 hypre_assert( hypre_ParVectorNumVectors(vector) == 1 );
848 if ( hypre_ParVectorNumVectors(vector) != 1 ) hypre_error_in_arg(1);
849
850 hypre_MPI_Comm_rank(comm, &myid);
851 hypre_MPI_Comm_size(comm, &num_procs);
852
853 hypre_sprintf(new_filename,"%s.%05d", filename, myid);
854
855 if ((file = fopen(new_filename, "w")) == NULL)
856 {
857 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Error: can't open output file %s\n");
858 return hypre_error_flag;
859 }
860
861 local_data = hypre_VectorData(hypre_ParVectorLocalVector(vector));
862
863 hypre_fprintf(file, "%b \n", global_size);
864 for (i=0; i < 2; i++)
865 {
866 hypre_fprintf(file, "%b ", partitioning[i] + base_j);
867 }
868 hypre_fprintf(file, "\n");
869
870 part0 = partitioning[0];
871 for (j = part0; j < partitioning[1]; j++)
872 {
873 hypre_fprintf(file, "%b %.14e\n", j + base_j, local_data[(HYPRE_Int)(j-part0)]);
874 }
875
876 fclose(file);
877
878 return hypre_error_flag;
879 }
880
881 /*--------------------------------------------------------------------------
882 * hypre_ParVectorReadIJ
883 * Warning: wrong base for assumed partition if base > 0
884 *--------------------------------------------------------------------------*/
885
886 HYPRE_Int
hypre_ParVectorReadIJ(MPI_Comm comm,const char * filename,HYPRE_Int * base_j_ptr,hypre_ParVector ** vector_ptr)887 hypre_ParVectorReadIJ( MPI_Comm comm,
888 const char *filename,
889 HYPRE_Int *base_j_ptr,
890 hypre_ParVector **vector_ptr )
891 {
892 HYPRE_BigInt global_size, J;
893 hypre_ParVector *vector;
894 hypre_Vector *local_vector;
895 HYPRE_Complex *local_data;
896 HYPRE_BigInt *partitioning;
897 HYPRE_Int base_j;
898
899 HYPRE_Int myid, num_procs, i, j;
900 char new_filename[255];
901 FILE *file;
902
903 hypre_MPI_Comm_size(comm, &num_procs);
904 hypre_MPI_Comm_rank(comm, &myid);
905
906 hypre_sprintf(new_filename,"%s.%05d", filename, myid);
907
908 if ((file = fopen(new_filename, "r")) == NULL)
909 {
910 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Error: can't open output file %s\n");
911 return hypre_error_flag;
912 }
913
914 hypre_fscanf(file, "%b", &global_size);
915 /* this may need to be changed so that the base is available in the file! */
916 partitioning = hypre_CTAlloc(HYPRE_BigInt, 2, HYPRE_MEMORY_HOST);
917
918 hypre_fscanf(file, "%b", partitioning);
919 for (i = 0; i < 2; i++)
920 {
921 hypre_fscanf(file, "%b", partitioning+i);
922 }
923 /* This is not yet implemented correctly! */
924 base_j = 0;
925 vector = hypre_ParVectorCreate(comm, global_size,
926 partitioning);
927
928 hypre_ParVectorInitialize(vector);
929
930 local_vector = hypre_ParVectorLocalVector(vector);
931 local_data = hypre_VectorData(local_vector);
932
933 for (j = 0; j < (HYPRE_Int)(partitioning[1] - partitioning[0]); j++)
934 {
935 hypre_fscanf(file, "%b %le", &J, local_data + j);
936 }
937
938 fclose(file);
939
940 *base_j_ptr = base_j;
941 *vector_ptr = vector;
942
943 /* multivector code not written yet */
944 hypre_assert( hypre_ParVectorNumVectors(vector) == 1 );
945 if ( hypre_ParVectorNumVectors(vector) != 1 ) hypre_error(HYPRE_ERROR_GENERIC);
946
947 return hypre_error_flag;
948 }
949
950 /*--------------------------------------------------------------------
951 * hypre_FillResponseParToVectorAll
952 * Fill response function for determining the send processors
953 * data exchange
954 *--------------------------------------------------------------------*/
955
956 HYPRE_Int
hypre_FillResponseParToVectorAll(void * p_recv_contact_buf,HYPRE_Int contact_size,HYPRE_Int contact_proc,void * ro,MPI_Comm comm,void ** p_send_response_buf,HYPRE_Int * response_message_size)957 hypre_FillResponseParToVectorAll( void *p_recv_contact_buf,
958 HYPRE_Int contact_size,
959 HYPRE_Int contact_proc,
960 void *ro,
961 MPI_Comm comm,
962 void **p_send_response_buf,
963 HYPRE_Int *response_message_size )
964 {
965 HYPRE_Int myid;
966 HYPRE_Int i, index, count, elength;
967
968 HYPRE_BigInt *recv_contact_buf = (HYPRE_BigInt * ) p_recv_contact_buf;
969
970 hypre_DataExchangeResponse *response_obj = (hypre_DataExchangeResponse*)ro;
971
972 hypre_ProcListElements *send_proc_obj = (hypre_ProcListElements*)response_obj->data2;
973 hypre_MPI_Comm_rank(comm, &myid );
974
975 /*check to see if we need to allocate more space in send_proc_obj for ids*/
976 if (send_proc_obj->length == send_proc_obj->storage_length)
977 {
978 send_proc_obj->storage_length +=10; /*add space for 10 more processors*/
979 send_proc_obj->id = hypre_TReAlloc(send_proc_obj->id, HYPRE_Int,
980 send_proc_obj->storage_length, HYPRE_MEMORY_HOST);
981 send_proc_obj->vec_starts =
982 hypre_TReAlloc(send_proc_obj->vec_starts, HYPRE_Int,
983 send_proc_obj->storage_length + 1, HYPRE_MEMORY_HOST);
984 }
985
986 /*initialize*/
987 count = send_proc_obj->length;
988 index = send_proc_obj->vec_starts[count]; /*this is the number of elements*/
989
990 /*send proc*/
991 send_proc_obj->id[count] = contact_proc;
992
993 /*do we need more storage for the elements?*/
994 if (send_proc_obj->element_storage_length < index + contact_size)
995 {
996 elength = hypre_max(contact_size, 10);
997 elength += index;
998 send_proc_obj->elements = hypre_TReAlloc(send_proc_obj->elements,
999 HYPRE_BigInt, elength, HYPRE_MEMORY_HOST);
1000 send_proc_obj->element_storage_length = elength;
1001 }
1002 /*populate send_proc_obj*/
1003 for (i=0; i< contact_size; i++)
1004 {
1005 send_proc_obj->elements[index++] = recv_contact_buf[i];
1006 }
1007 send_proc_obj->vec_starts[count+1] = index;
1008 send_proc_obj->length++;
1009
1010 /*output - no message to return (confirmation) */
1011 *response_message_size = 0;
1012
1013 return hypre_error_flag;
1014 }
1015
1016 /* -----------------------------------------------------------------------------
1017 * return the sum of all local elements of the vector
1018 * ----------------------------------------------------------------------------- */
1019
hypre_ParVectorLocalSumElts(hypre_ParVector * vector)1020 HYPRE_Complex hypre_ParVectorLocalSumElts( hypre_ParVector * vector )
1021 {
1022 return hypre_SeqVectorSumElts( hypre_ParVectorLocalVector(vector) );
1023 }
1024
1025 HYPRE_Int
hypre_ParVectorGetValuesHost(hypre_ParVector * vector,HYPRE_Int num_values,HYPRE_BigInt * indices,HYPRE_BigInt base,HYPRE_Complex * values)1026 hypre_ParVectorGetValuesHost(hypre_ParVector *vector,
1027 HYPRE_Int num_values,
1028 HYPRE_BigInt *indices,
1029 HYPRE_BigInt base,
1030 HYPRE_Complex *values)
1031 {
1032 HYPRE_Int i, ierr = 0;
1033 HYPRE_BigInt first_index = hypre_ParVectorFirstIndex(vector);
1034 HYPRE_BigInt last_index = hypre_ParVectorLastIndex(vector);
1035 hypre_Vector *local_vector = hypre_ParVectorLocalVector(vector);
1036 HYPRE_Complex *data = hypre_VectorData(local_vector);
1037
1038 /*
1039 if (hypre_VectorOwnsData(local_vector) == 0)
1040 {
1041 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Vector does not own data! -- hypre_ParVectorGetValues.");
1042 return hypre_error_flag;
1043 }
1044 */
1045
1046 if (indices)
1047 {
1048 #ifdef HYPRE_USING_OPENMP
1049 #pragma omp parallel for private(i) reduction(+:ierr) HYPRE_SMP_SCHEDULE
1050 #endif
1051 for (i = 0; i < num_values; i++)
1052 {
1053 HYPRE_BigInt index = indices[i] - base;
1054 if (index < first_index || index > last_index)
1055 {
1056 ierr ++;
1057 }
1058 else
1059 {
1060 HYPRE_Int local_index = (HYPRE_Int) (index - first_index);
1061 values[i] = data[local_index];
1062 }
1063 }
1064
1065 if (ierr)
1066 {
1067 hypre_error_in_arg(3);
1068 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Index out of range! -- hypre_ParVectorGetValues.");
1069 hypre_printf("Index out of range! -- hypre_ParVectorGetValues\n");
1070 }
1071 }
1072 else
1073 {
1074 if (num_values > hypre_VectorSize(local_vector))
1075 {
1076 hypre_error_in_arg(2);
1077 return hypre_error_flag;
1078 }
1079
1080 #ifdef HYPRE_USING_OPENMP
1081 #pragma omp parallel for private(i) HYPRE_SMP_SCHEDULE
1082 #endif
1083 for (i = 0; i < num_values; i++)
1084 {
1085 values[i] = data[i];
1086 }
1087 }
1088
1089 return hypre_error_flag;
1090 }
1091
1092 HYPRE_Int
hypre_ParVectorGetValues2(hypre_ParVector * vector,HYPRE_Int num_values,HYPRE_BigInt * indices,HYPRE_BigInt base,HYPRE_Complex * values)1093 hypre_ParVectorGetValues2(hypre_ParVector *vector,
1094 HYPRE_Int num_values,
1095 HYPRE_BigInt *indices,
1096 HYPRE_BigInt base,
1097 HYPRE_Complex *values)
1098 {
1099 #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
1100 if (HYPRE_EXEC_DEVICE == hypre_GetExecPolicy1( hypre_ParVectorMemoryLocation(vector) ))
1101 {
1102 hypre_ParVectorGetValuesDevice(vector, num_values, indices, base, values);
1103 }
1104 else
1105 #endif
1106 {
1107 hypre_ParVectorGetValuesHost(vector, num_values, indices, base, values);
1108 }
1109
1110 return hypre_error_flag;
1111 }
1112
1113 HYPRE_Int
hypre_ParVectorGetValues(hypre_ParVector * vector,HYPRE_Int num_values,HYPRE_BigInt * indices,HYPRE_Complex * values)1114 hypre_ParVectorGetValues(hypre_ParVector *vector,
1115 HYPRE_Int num_values,
1116 HYPRE_BigInt *indices,
1117 HYPRE_Complex *values)
1118 {
1119 return hypre_ParVectorGetValues2(vector, num_values, indices, 0, values);
1120 }
1121