/[escript]/branches/symbolic_from_3470/paso/src/AMG_Prolongation.c
ViewVC logotype

Diff of /branches/symbolic_from_3470/paso/src/AMG_Prolongation.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3867 by caltinay, Tue Jan 31 04:55:05 2012 UTC revision 3868 by caltinay, Thu Mar 15 06:07:08 2012 UTC
# Line 1  Line 1 
   
1  /*******************************************************  /*******************************************************
2  *  *
3  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
# Line 47  Line 46 
46        
47  */  */
48    
49    Paso_SystemMatrix* Paso_Preconditioner_AMG_getProlongation(Paso_SystemMatrix* A_p,
 Paso_SparseMatrix* Paso_Preconditioner_AMG_getProlongation(Paso_SparseMatrix* A_p,  
50                                                             const index_t* offset_S, const dim_t* degree_S, const index_t* S,                                                             const index_t* offset_S, const dim_t* degree_S, const index_t* S,
51                                 const dim_t n_C, const index_t* counter_C, const index_t interpolation_method)                                 const dim_t n_C, index_t* counter_C, const index_t interpolation_method)
52  {  {
53     Paso_SparseMatrix* out=NULL;     Esys_MPIInfo *mpi_info=Esys_MPIInfo_getReference(A_p->mpi_info);
54     Paso_Pattern *outpattern=NULL;     Paso_SystemMatrix *out=NULL;
55     const dim_t n_block=A_p->row_block_size;     Paso_SystemMatrixPattern *pattern=NULL;
56     index_t *ptr=NULL, *index=NULL,j, iptr;     Paso_Distribution *input_dist=NULL, *output_dist=NULL;
57     const dim_t n =A_p->numRows;     Paso_SharedComponents *send =NULL, *recv=NULL;
58     dim_t i,p,z, len_P;     Paso_Connector *col_connector=NULL;
59         Paso_Pattern *main_pattern=NULL, *couple_pattern=NULL;
60     ptr=MEMALLOC(n+1,index_t);     const dim_t row_block_size=A_p->row_block_size;
61     if (! Esys_checkPtr(ptr)) {     const dim_t col_block_size=A_p->col_block_size;
62       const dim_t my_n=A_p->mainBlock->numCols;
63       const dim_t overlap_n=A_p->col_coupleBlock->numCols;
64       const dim_t num_threads=omp_get_max_threads();
65       index_t size=mpi_info->size, *dist=NULL;
66       index_t *main_p=NULL, *couple_p=NULL, *main_idx=NULL, *couple_idx=NULL;
67       index_t *shared=NULL, *offsetInShared=NULL;
68       index_t *recv_shared=NULL, *send_shared=NULL;
69       index_t sum, i, j, k, l, p, q, iptr;
70       index_t my_n_C, global_label, num_neighbors;
71       #ifdef ESYS_MPI
72       index_t rank=mpi_info->rank;
73       #endif
74       Esys_MPI_rank *neighbor=NULL;
75       #ifdef ESYS_MPI
76         MPI_Request* mpi_requests=NULL;
77         MPI_Status* mpi_stati=NULL;
78       #else
79         int *mpi_requests=NULL, *mpi_stati=NULL;
80       #endif
81    
82       /* number of C points in current distribution */
83       my_n_C = 0;
84       sum=0;
85       if (num_threads>1) {
86         #pragma omp parallel private(i,sum)
87         {
88        sum=0;
89        #pragma omp for schedule(static)
90        for (i=0;i<my_n;++i) {
91          if (counter_C[i] != -1) {
92            sum++;
93          }
94        }
95        #pragma omp critical
96        {
97            my_n_C += sum;
98        }
99         }
100       } else { /* num_threads=1 */
101         for (i=0;i<my_n;++i) {
102             if (counter_C[i] != -1) {
103                my_n_C++;
104             }
105          }
106       }
107    
108             /* create row distribution (output_distribution) and col distribution
109        /* count the number of entries per row in the Prolongation matrix :*/        (input_distribution) */
110         /* ??? should I alloc an new Esys_MPIInfo object or reuse the one in
111        #pragma omp parallel for private(i,z,iptr,j,p) schedule(static)        system matrix A. for now, I'm reuse A->mpi_info ??? */
112        for (i=0;i<n;++i) {     dist = A_p->pattern->output_distribution->first_component;
113       if (counter_C[i]>=0) {     output_dist=Paso_Distribution_alloc(mpi_info, dist, 1, 0);
114          z=1;    /* i is a C unknown */     dist = TMPMEMALLOC(size+1, index_t); /* now prepare for col distribution */
115       } else {     Esys_checkPtr(dist);
116          z=0;     #ifdef ESYS_MPI
117          iptr=offset_S[i];     MPI_Allgather(&my_n_C, 1, MPI_INT, dist, 1, MPI_INT, mpi_info->comm);
118          for (p=0; p<degree_S[i]; ++p) {     #endif
119             j=S[iptr+p];  /* this is a strong connection */     global_label=0;
120             if (counter_C[j]>=0) z++; /* and is in C */     for (i=0; i<size; i++) {
121         k = dist[i];
122         dist[i] = global_label;
123         global_label += k;
124       }
125       dist[size] = global_label;
126    
127       input_dist=Paso_Distribution_alloc(mpi_info, dist, 1, 0);
128       TMPMEMFREE(dist);
129    
130       /* create pattern for mainBlock and coupleBlock */
131       main_p = MEMALLOC(my_n+1, index_t);
132       couple_p = MEMALLOC(my_n+1, index_t);
133       if (!(Esys_checkPtr(main_p) || Esys_checkPtr(couple_p))) {
134         /* count the number of entries per row in the Prolongation matrix :*/
135         #pragma omp parallel for private(i,l,k,iptr,j,p) schedule(static)
136         for (i=0; i<my_n; i++) {
137        l = 0;
138        if (counter_C[i]>=0) {
139          k = 1;    /* i is a C unknown */
140        } else {
141          k = 0;
142          iptr = offset_S[i];
143          for (p=0; p<degree_S[i]; p++) {
144            j = S[iptr+p];  /* this is a strong connection */
145            if (counter_C[j]>=0) { /* and is in C */
146            if (j <my_n) k++;
147            else {
148              l++;
149            }
150          }          }
151       }        }
152       ptr[i]=z;      }
153        }      main_p[i] = k;
154        len_P=Paso_Util_cumsum(n,ptr);      couple_p[i] = l;
155        ptr[n]=len_P;       }
156          
157        /* allocate and create index vector for prolongation: */       /* number of unknowns in the col-coupleBlock of the interplation matrix */
158        index=MEMALLOC(len_P,index_t);       sum = 0;
159           for (i=0;i<overlap_n;++i) {
160        if (! Esys_checkPtr(index)) {      if (counter_C[i+my_n] > -1) {
161       #pragma omp parallel for private(i,z,iptr,j,p)  schedule(static)        counter_C[i+my_n] -= my_n_C;
162       for (i=0;i<n;++i) {        sum++;
163          if (counter_C[i]>=0) {      }
164             index[ptr[i]]=counter_C[i];  /* i is a C unknown */       }
165          } else {  
166             z=0;       /* allocate and create index vector for prolongation: */
167             iptr=offset_S[i];       p = Paso_Util_cumsum(my_n, main_p);
168             for (p=0; p<degree_S[i]; ++p) {       main_p[my_n] = p;
169            j=S[iptr+p];  /* this is a strong connection */       main_idx = MEMALLOC(p, index_t);
170            if (counter_C[j]>=0) {  /* and is in C */       p = Paso_Util_cumsum(my_n, couple_p);
171               index[ptr[i]+z]=counter_C[j];       couple_p[my_n] = p;
172               z++; /* and is in C */       couple_idx = MEMALLOC(p, index_t);
173            }       if (!(Esys_checkPtr(main_idx) || Esys_checkPtr(couple_idx))) {
174             }      #pragma omp parallel for private(i,k,l,iptr,j,p)  schedule(static)
175        for (i=0; i<my_n; i++) {
176          if (counter_C[i]>=0) {
177            main_idx[main_p[i]]=counter_C[i];  /* i is a C unknown */
178          } else {
179            k = 0;
180            l = 0;
181            iptr = offset_S[i];
182            for (p=0; p<degree_S[i]; p++) {
183              j = S[iptr+p]; /* this is a strong connection */
184              if (counter_C[j] >=0) { /* and is in C */
185            if (j < my_n) {
186              main_idx[main_p[i]+k] = counter_C[j];
187              k++;
188            } else {
189              couple_idx[couple_p[i]+l] = counter_C[j];
190              l++;
191            }
192              }
193          }          }
194       }        }
195        }      }
196     }         }
197     if (Esys_noError()) {     }
198       outpattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,n,n_C,ptr,index);  
199       if (Esys_noError()) {  
200         main_pattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, my_n,
201                my_n_C, main_p, main_idx);
202         couple_pattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, my_n,
203                sum, couple_p, couple_idx);
204     } else {     } else {
205        MEMFREE(ptr);       MEMFREE(main_p);
206        MEMFREE(index);       MEMFREE(main_idx);
207         MEMFREE(couple_p);
208         MEMFREE(couple_idx);
209       }
210    
211       /* prepare the receiver for the col_connector.
212          Note that the allocation for "shared" assumes the send and receive buffer
213          of the interpolation matrix P is no larger than that of matrix A_p. */
214       neighbor = TMPMEMALLOC(size, Esys_MPI_rank);
215       offsetInShared = TMPMEMALLOC(size+1, index_t);
216       recv = A_p->col_coupler->connector->recv;
217       send = A_p->col_coupler->connector->send;
218       i = recv->numSharedComponents;
219       recv_shared = TMPMEMALLOC(i,index_t);
220       memset(recv_shared, 0, sizeof(index_t)*i);
221       k = send->numSharedComponents;
222       send_shared = TMPMEMALLOC(k,index_t);
223       if (k > i) i = k;
224       shared = TMPMEMALLOC(i, index_t);
225    
226       #ifdef ESYS_MPI
227         mpi_requests=TMPMEMALLOC(size*2,MPI_Request);
228         mpi_stati=TMPMEMALLOC(size*2,MPI_Status);
229       #else
230         mpi_requests=TMPMEMALLOC(size*2,int);
231         mpi_stati=TMPMEMALLOC(size*2,int);
232       #endif
233    
234       for (p=0; p<send->numNeighbors; p++) {
235         i = send->offsetInShared[p];
236         #ifdef ESYS_MPI
237         MPI_Irecv (&(send_shared[i]), send->offsetInShared[p+1]-i, MPI_INT,
238            send->neighbor[p], mpi_info->msg_tag_counter+send->neighbor[p],
239            mpi_info->comm, &mpi_requests[p]);
240         #endif
241       }
242    
243       num_neighbors = 0;
244       q = 0;
245       p = recv->numNeighbors;
246       offsetInShared[0]=0;
247       for (i=0; i<p; i++) {
248         l = 0;
249         k = recv->offsetInShared[i+1];
250         for (j=recv->offsetInShared[i]; j<k; j++) {
251        if (counter_C[recv->shared[j]] > -1) {
252          shared[q] = my_n_C + q;
253          recv_shared[recv->shared[j]-my_n] = 1;
254          q++;
255          l = 1;
256        }
257         }
258         if (l == 1) {
259        iptr = recv->neighbor[i];
260        neighbor[num_neighbors] = iptr;
261        num_neighbors++;
262        offsetInShared[num_neighbors] = q;
263         }
264         #ifdef ESYS_MPI
265         MPI_Issend(&(recv_shared[recv->offsetInShared[i]]),
266            k-recv->offsetInShared[i], MPI_INT, recv->neighbor[i],
267            mpi_info->msg_tag_counter+rank, mpi_info->comm,
268            &mpi_requests[i+send->numNeighbors]);
269         #endif
270       }
271       recv = Paso_SharedComponents_alloc(my_n_C, num_neighbors, neighbor, shared,
272                                          offsetInShared, 1, 0, mpi_info);
273    
274       /* now we can build the sender */
275       #ifdef ESYS_MPI
276       MPI_Waitall(recv->numNeighbors+send->numNeighbors, mpi_requests, mpi_stati);
277       #endif
278       mpi_info->msg_tag_counter += size;
279       TMPMEMFREE(mpi_requests);
280       TMPMEMFREE(mpi_stati);
281    
282       num_neighbors = 0;
283       q = 0;
284       p = send->numNeighbors;
285       offsetInShared[0]=0;
286       for (i=0; i<p; i++) {
287         l = 0;
288         k = send->offsetInShared[i+1];
289         for (j=send->offsetInShared[i]; j<k; j++) {
290        if (send_shared[j] == 1) {
291              shared[q] = counter_C[send->shared[j]];  
292              q++;
293              l = 1;
294            }
295         }
296         if (l == 1) {
297            iptr = send->neighbor[i];
298            neighbor[num_neighbors] = iptr;
299            num_neighbors++;
300            offsetInShared[num_neighbors] = q;
301         }
302     }     }
303     /* now we need to create a matrix and fill it */  
304       send = Paso_SharedComponents_alloc(my_n_C, num_neighbors, neighbor, shared,
305                          offsetInShared, 1, 0, mpi_info);
306       col_connector = Paso_Connector_alloc(send, recv);
307       Paso_SharedComponents_free(recv);
308       Paso_SharedComponents_free(send);
309       TMPMEMFREE(recv_shared);
310       TMPMEMFREE(send_shared);
311       TMPMEMFREE(neighbor);
312       TMPMEMFREE(offsetInShared);
313       TMPMEMFREE(shared);
314    
315       /* now we need to create the System Matrix
316          TO BE FIXED: at this stage, we only construction col_couple_pattern
317          and col_connector for interpolation matrix P. To be completed,
318          row_couple_pattern and row_connector need to be constructed as well */
319     if (Esys_noError()) {     if (Esys_noError()) {
320        out=Paso_SparseMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK,outpattern,n_block,n_block,FALSE);       pattern = Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,
321            output_dist, input_dist, main_pattern, couple_pattern,
322            couple_pattern, col_connector, col_connector);
323         out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DIAGONAL_BLOCK, pattern,
324            row_block_size, col_block_size, FALSE);
325     }     }
326      
327       /* now fill in the matrix */
328     if (Esys_noError()) {     if (Esys_noError()) {
329        if ( (interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING) || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {       if ((interpolation_method == PASO_CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
330          if (n_block == 1) {          || ( interpolation_method == PASO_CLASSIC_INTERPOLATION) ) {
331          Paso_Preconditioner_AMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);          if (row_block_size == 1) {
332          } else {            Paso_Preconditioner_AMG_setClassicProlongation(out, A_p, offset_S, degree_S, S, counter_C);
333          Paso_Preconditioner_AMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);          } else {
334          }            Paso_Preconditioner_AMG_setClassicProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
335        } else {          }
336          if (n_block == 1) {       } else {
337              Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, counter_C);          if (row_block_size == 1) {
338          } else {            Paso_Preconditioner_AMG_setDirectProlongation(out, A_p, offset_S, degree_S, S, counter_C);
339          Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, counter_C);          } else {
340          }            Paso_Preconditioner_AMG_setDirectProlongation_Block(out, A_p, offset_S, degree_S, S, counter_C);
341        }          }
342     }       }
343     Paso_Pattern_free(outpattern);     }  
344    
345       /* clean up */
346       Paso_SystemMatrixPattern_free(pattern);
347       Paso_Pattern_free(main_pattern);
348       Paso_Pattern_free(couple_pattern);
349       Paso_Connector_free(col_connector);
350       Paso_Distribution_free(output_dist);
351       Paso_Distribution_free(input_dist);
352     if (Esys_noError()) {     if (Esys_noError()) {
353        return out;        return out;
354     } else {     } else {
355        Paso_SparseMatrix_free(out);        Paso_SystemMatrix_free(out);
356        return NULL;        return NULL;
357     }     }
358        
# Line 155  Paso_SparseMatrix* Paso_Preconditioner_A Line 373  Paso_SparseMatrix* Paso_Preconditioner_A
373    
374  */  */
375    
376  void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SparseMatrix* P_p,  void Paso_Preconditioner_AMG_setDirectProlongation(Paso_SystemMatrix* P,
377                                 const Paso_SparseMatrix* A_p,      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
378                             const index_t *counter_C) {          const index_t* S, const index_t *counter_C) {
379       Paso_SparseMatrix *main_block=P->mainBlock;
380       Paso_SparseMatrix *couple_block=P->col_coupleBlock;
381       Paso_Pattern *main_pattern=main_block->pattern;
382       Paso_Pattern *couple_pattern=couple_block->pattern;
383       const dim_t my_n=A->mainBlock->numRows;
384       index_t range;
385    
386     dim_t i;     dim_t i;
    const dim_t n =A_p->numRows;  
387     register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;     register double alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos, A_ij, A_ii, rtmp;
388     register index_t iPtr, j, offset;     register index_t iPtr, j, offset;
389     index_t *where_p, *start_p;     index_t *where_p, *start_p;
390        
391     #pragma omp parallel for private(A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij , rtmp)  schedule(static)     #pragma omp parallel for private(i,offset,sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos,A_ii,range,iPtr,j,A_ij,start_p,where_p,alpha,beta,rtmp) schedule(static)
392     for (i=0;i<n;++i) {     for (i=0; i<my_n; i++) {
393        if (counter_C[i]>=0) {        if (counter_C[i]>=0) {
394          offset = P_p->pattern->ptr[i];          offset = main_pattern->ptr[i];
395          P_p->val[offset]=1.;  /* i is a C row */          main_block->val[offset]=1.;  /* i is a C row */
396        } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {        } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
397             (couple_pattern->ptr[i +1] > couple_pattern->ptr[i])) {
398       /* if i is an F row we first calculate alpha and beta: */       /* if i is an F row we first calculate alpha and beta: */
399    
400       sum_all_neg=0; /* sum of all negative values in row i of A */       sum_all_neg=0; /* sum of all negative values in row i of A */
401       sum_all_pos=0; /* sum of all positive values in row i of A */       sum_all_pos=0; /* sum of all positive values in row i of A */
402       sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/       sum_strong_neg=0; /* sum of all negative values A_ij where j is in C and strongly connected to i*/
403       sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/       sum_strong_pos=0; /* sum of all positive values A_ij where j is in C and strongly connected to i*/
404       A_ii=0;       A_ii=0;
405       for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {  
406          j=A_p->pattern->index[iPtr];       /* first check the mainBlock */
407          A_ij=A_p->val[iPtr];       range = A->mainBlock->pattern->ptr[i + 1];
408         for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
409            j=A->mainBlock->pattern->index[iPtr];
410            A_ij=A->mainBlock->val[iPtr];
411          if(j==i) {          if(j==i) {
412             A_ii=A_ij;             A_ii=A_ij;
413          } else {          } else {
# Line 191  void Paso_Preconditioner_AMG_setDirectPr Line 420  void Paso_Preconditioner_AMG_setDirectPr
420                        
421             if (counter_C[j]>=0) {             if (counter_C[j]>=0) {
422            /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */            /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
423            start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);            start_p=&(main_pattern->index[main_pattern->ptr[i]]);
424            where_p=(index_t*)bsearch(&(counter_C[j]), start_p,            where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
425                          P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],                          main_pattern->ptr[i + 1] - main_pattern->ptr[i],
426                          sizeof(index_t),                          sizeof(index_t),
427                          Paso_comparIndex);                          Paso_comparIndex);
428            if (! (where_p == NULL) ) { /* yes i strongly connected with j */            if (! (where_p == NULL) ) { /* yes i strongly connected with j */
429              offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);              offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p);
430              P_p->val[offset]=A_ij; /* will be modified later */              main_block->val[offset]=A_ij; /* will be modified later */
431              if (A_ij< 0)  {              if (A_ij< 0)  {
432                 sum_strong_neg+=A_ij;                 sum_strong_neg+=A_ij;
433              } else {              } else {
# Line 209  void Paso_Preconditioner_AMG_setDirectPr Line 438  void Paso_Preconditioner_AMG_setDirectPr
438    
439          }          }
440       }       }
441    
442             /* now we deal with the col_coupleBlock */
443             range = A->col_coupleBlock->pattern->ptr[i + 1];
444             for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
445                j=A->col_coupleBlock->pattern->index[iPtr];
446                A_ij=A->col_coupleBlock->val[iPtr];
447                if(A_ij< 0)  {
448                      sum_all_neg+=A_ij;
449                } else {
450                      sum_all_pos+=A_ij;
451                }
452    
453                if (counter_C[j+my_n]>=0) {
454                      /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
455                      start_p=&(couple_pattern->index[couple_pattern->ptr[i]]);
456                      where_p=(index_t*)bsearch(&(counter_C[j+my_n]), start_p,
457                                                couple_pattern->ptr[i + 1] - couple_pattern->ptr[i],
458                                                sizeof(index_t),
459                                                Paso_comparIndex);
460                      if (! (where_p == NULL) ) { /* yes i stronly connect with j */
461                            offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p);
462                            couple_block->val[offset]=A_ij; /* will be modified later */
463                            if (A_ij< 0)  {
464                               sum_strong_neg+=A_ij;
465                            } else {
466                               sum_strong_pos+=A_ij;
467                            }
468                      }
469                }
470             }
471    
472       if(sum_strong_neg<0) {       if(sum_strong_neg<0) {
473          alpha= sum_all_neg/sum_strong_neg;          alpha= sum_all_neg/sum_strong_neg;
474       } else {       } else {
# Line 225  void Paso_Preconditioner_AMG_setDirectPr Line 485  void Paso_Preconditioner_AMG_setDirectPr
485          alpha*=rtmp;          alpha*=rtmp;
486          beta*=rtmp;          beta*=rtmp;
487       }       }
488       for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {  
489          A_ij=P_p->val[iPtr];       range = main_pattern->ptr[i + 1];
490         for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
491            A_ij=main_block->val[iPtr];
492          if (A_ij > 0 ) {          if (A_ij > 0 ) {
493             P_p->val[iPtr]*=beta;             main_block->val[iPtr] = A_ij * beta;
494          } else {          } else {
495             P_p->val[iPtr]*=alpha;             main_block->val[iPtr] = A_ij * alpha;
496          }          }
497       }       }
498    
499             range = couple_pattern->ptr[i + 1];
500             for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
501                A_ij=couple_block->val[iPtr];
502                if (A_ij > 0 ) {
503                   couple_block->val[iPtr] = A_ij * beta;
504                } else {
505                   couple_block->val[iPtr] = A_ij * alpha;
506                }
507             }
508        }        }
509     }     }
510  }  }
511    
512  void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SparseMatrix* P_p,  void Paso_Preconditioner_AMG_setDirectProlongation_Block(Paso_SystemMatrix* P,
513                                   const Paso_SparseMatrix* A_p,          Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
514                                   const index_t *counter_C) {          const index_t* S, const index_t *counter_C) {
515       Paso_SparseMatrix *main_block=P->mainBlock;
516       Paso_SparseMatrix *couple_block=P->col_coupleBlock;
517       Paso_Pattern *main_pattern=main_block->pattern;
518       Paso_Pattern *couple_pattern=couple_block->pattern;
519       const dim_t row_block_size=A->row_block_size;
520       const dim_t A_block = A->block_size;
521       const dim_t my_n=A->mainBlock->numRows;
522       index_t range;
523    
524     dim_t i;     dim_t i;
    const dim_t n =A_p->numRows;  
    const dim_t row_block=A_p->row_block_size;  
    const dim_t A_block = A_p->block_size;  
525     double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;     double *alpha, *beta, *sum_all_neg, *sum_all_pos, *sum_strong_neg, *sum_strong_pos, *A_ii;
526     register double A_ij, rtmp;     register double A_ij, rtmp;
527     register index_t iPtr, j, offset, ib;     register index_t iPtr, j, offset, ib;
528     index_t *where_p, *start_p;     index_t *where_p, *start_p;
529        
530     #pragma omp parallel private(ib, rtmp, A_ii, offset, where_p, start_p, i, alpha, beta, sum_all_neg, sum_all_pos, sum_strong_neg, sum_strong_pos,iPtr,j, A_ij )       #pragma omp parallel private(i,offset,ib,sum_all_neg,sum_all_pos,sum_strong_neg,sum_strong_pos,A_ii,range,iPtr,j,A_ij,start_p,where_p,alpha,beta,rtmp)
531     {     {
532        sum_all_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values in row i of A */        sum_all_neg=TMPMEMALLOC(row_block_size, double); /* sum of all negative values in row i of A */
533        sum_all_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values in row i of A */        sum_all_pos=TMPMEMALLOC(row_block_size, double); /* sum of all positive values in row i of A */
534        sum_strong_neg=TMPMEMALLOC(row_block, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/        sum_strong_neg=TMPMEMALLOC(row_block_size, double); /* sum of all negative values A_ij where j is in C and strongly connected to i*/
535        sum_strong_pos=TMPMEMALLOC(row_block, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/        sum_strong_pos=TMPMEMALLOC(row_block_size, double); /* sum of all positive values A_ij where j is in C and strongly connected to i*/
536        alpha=TMPMEMALLOC(row_block, double);        alpha=TMPMEMALLOC(row_block_size, double);
537        beta=TMPMEMALLOC(row_block, double);        beta=TMPMEMALLOC(row_block_size, double);
538        A_ii=TMPMEMALLOC(row_block, double);        A_ii=TMPMEMALLOC(row_block_size, double);
539                
540        #pragma omp for schedule(static)        #pragma omp for schedule(static)
541        for (i=0;i<n;++i) {        for (i=0;i<my_n;++i) {
542       if (counter_C[i]>=0) {       if (counter_C[i]>=0) { /* i is a C row */
543          offset = P_p->pattern->ptr[i];          offset = main_pattern->ptr[i];
544          for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.;  /* i is a C row */          for (ib =0; ib<row_block_size; ++ib)
545           } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {          main_block->val[row_block_size*offset+ib]=1.;
546             } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
547                (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
548          /* if i is an F row we first calculate alpha and beta: */          /* if i is an F row we first calculate alpha and beta: */
549          for (ib =0; ib<row_block; ++ib) {          for (ib =0; ib<row_block_size; ++ib) {
550             sum_all_neg[ib]=0;             sum_all_neg[ib]=0;
551             sum_all_pos[ib]=0;             sum_all_pos[ib]=0;
552             sum_strong_neg[ib]=0;             sum_strong_neg[ib]=0;
553             sum_strong_pos[ib]=0;             sum_strong_pos[ib]=0;
554             A_ii[ib]=0;             A_ii[ib]=0;
555          }          }
556          for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {  
557             j=A_p->pattern->index[iPtr];          /* first check the mainBlock */
558            range = A->mainBlock->pattern->ptr[i + 1];
559            for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
560               j=A->mainBlock->pattern->index[iPtr];
561             if(j==i) {             if(j==i) {
562            for (ib =0; ib<row_block; ++ib) A_ii[ib]=A_p->val[A_block*iPtr+ib+row_block*ib];            for (ib =0; ib<row_block_size; ++ib)
563                A_ii[ib]=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
564             } else {             } else {
565            for (ib =0; ib<row_block; ++ib) {            for (ib =0; ib<row_block_size; ++ib) {
566               A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];               A_ij=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
567               if(A_ij< 0)  {               if(A_ij< 0)  {
568              sum_all_neg[ib]+=A_ij;              sum_all_neg[ib]+=A_ij;
569               } else {               } else {
# Line 289  void Paso_Preconditioner_AMG_setDirectPr Line 573  void Paso_Preconditioner_AMG_setDirectPr
573                        
574            if (counter_C[j]>=0) {            if (counter_C[j]>=0) {
575               /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */               /* is i strongly connected with j? We search for counter_C[j] in P[i,:] */
576               start_p=&(P_p->pattern->index[P_p->pattern->ptr[i]]);               start_p=&(main_pattern->index[main_pattern->ptr[i]]);
577               where_p=(index_t*)bsearch(&(counter_C[j]), start_p,               where_p=(index_t*)bsearch(&(counter_C[j]), start_p,
578                           P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i],                           main_pattern->ptr[i + 1]-main_pattern->ptr[i],
579                           sizeof(index_t),                           sizeof(index_t),
580                           Paso_comparIndex);                           Paso_comparIndex);
581               if (! (where_p == NULL) ) { /* yes i strongly connected with j */               if (! (where_p == NULL) ) { /* yes i strongly connected with j */
582                    offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                    offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p);
583                    for (ib =0; ib<row_block; ++ib) {                    for (ib =0; ib<row_block_size; ++ib) {
584                   A_ij=A_p->val[A_block*iPtr+ib+row_block*ib];                   A_ij=A->mainBlock->val[A_block*iPtr+ib+row_block_size*ib];
585                   P_p->val[row_block*offset+ib]=A_ij; /* will be modified later */                   main_block->val[row_block_size*offset+ib]=A_ij; /* will be modified later */
586                       if (A_ij< 0)  {                       if (A_ij< 0)  {
587                       sum_strong_neg[ib]+=A_ij;                       sum_strong_neg[ib]+=A_ij;
588                   } else {                   } else {
# Line 310  void Paso_Preconditioner_AMG_setDirectPr Line 594  void Paso_Preconditioner_AMG_setDirectPr
594                    
595             }             }
596          }          }
597          for (ib =0; ib<row_block; ++ib) {  
598            /* now we deal with the col_coupleBlock */
599            range = A->col_coupleBlock->pattern->ptr[i + 1];
600                for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
601                      j=A->col_coupleBlock->pattern->index[iPtr];
602                      for (ib =0; ib<row_block_size; ++ib) {
603                         A_ij=A->col_coupleBlock->val[A_block*iPtr+ib+row_block_size*ib];
604                         if(A_ij< 0)  {
605                            sum_all_neg[ib]+=A_ij;
606                         } else {
607                            sum_all_pos[ib]+=A_ij;
608                         }
609                      }
610    
611                      if (counter_C[j+my_n]>=0) {
612                         /* is i stronly connect with j? We serach for counter_C[j] in P[i,:] */
613                         start_p=&(couple_pattern->index[couple_pattern->ptr[i]]);
614                         where_p=(index_t*)bsearch(&(counter_C[j+my_n]), start_p,
615                                                 couple_pattern->ptr[i + 1]-couple_pattern->ptr[i],
616                                                 sizeof(index_t),
617                                                 Paso_comparIndex);
618                         if (! (where_p == NULL) ) { /* yes i stronly connect with j */
619                                  offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p);
620                                  for (ib =0; ib<row_block_size; ++ib) {
621                                     A_ij=A->col_coupleBlock->val[A_block*iPtr+ib+row_block_size*ib];
622                                     couple_block->val[row_block_size*offset+ib]=A_ij; /* will be modified later */
623    
624                                     if (A_ij< 0)  {
625                                         sum_strong_neg[ib]+=A_ij;
626                                     } else {
627                                        sum_strong_pos[ib]+=A_ij;
628                                     }
629                                  }
630                         }
631                      }
632                }
633    
634            for (ib =0; ib<row_block_size; ++ib) {
635             if(sum_strong_neg[ib]<0) {             if(sum_strong_neg[ib]<0) {
636            alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];            alpha[ib]= sum_all_neg[ib]/sum_strong_neg[ib];
637             } else {             } else {
# Line 328  void Paso_Preconditioner_AMG_setDirectPr Line 649  void Paso_Preconditioner_AMG_setDirectPr
649            beta[ib]*=rtmp;            beta[ib]*=rtmp;
650             }             }
651          }          }
652          
653          for (iPtr=P_p->pattern->ptr[i];iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {          range = main_pattern->ptr[i + 1];      
654             for (ib =0; ib<row_block; ++ib) {          for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
655            A_ij=P_p->val[row_block*iPtr+ib];             for (ib =0; ib<row_block_size; ++ib) {
656              A_ij=main_block->val[row_block_size*iPtr+ib];
657            if (A_ij > 0 ) {            if (A_ij > 0 ) {
658               P_p->val[row_block*iPtr+ib]*=beta[ib];               main_block->val[row_block_size*iPtr+ib] = A_ij * beta[ib];
659            } else {            } else {
660               P_p->val[row_block*iPtr+ib]*=alpha[ib];               main_block->val[row_block_size*iPtr+ib] = A_ij * alpha[ib];
661            }            }
662             }             }
663          }          }
664    
665            range = couple_pattern->ptr[i + 1];
666            for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
667              for (ib =0; ib<row_block_size; ++ib) {
668                      A_ij=couple_block->val[row_block_size*iPtr+ib];
669                      if (A_ij > 0 ) {
670                         couple_block->val[row_block_size*iPtr+ib] = A_ij * beta[ib];
671                      } else {
672                         couple_block->val[row_block_size*iPtr+ib] = A_ij * alpha[ib];
673                      }
674                   }
675                }
676    
677    
678       }       }
679        }/* end i loop */        }/* end i loop */
680        TMPMEMFREE(sum_all_neg);        TMPMEMFREE(sum_all_neg);
# Line 367  void Paso_Preconditioner_AMG_setDirectPr Line 703  void Paso_Preconditioner_AMG_setDirectPr
703                                
704    
705  */  */
706  void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SparseMatrix* P_p,  void Paso_Preconditioner_AMG_setClassicProlongation(Paso_SystemMatrix* P,
707                                 Paso_SparseMatrix* A_p,      Paso_SystemMatrix* A, const index_t* offset_S, const dim_t* degree_S,
708                                                     const index_t* offset_S, const dim_t* degree_S, const index_t* S,      const index_t* S, const index_t *counter_C) {
709                             const index_t *counter_C) {     Paso_SparseMatrix *main_block=P->mainBlock;
710       Paso_SparseMatrix *couple_block=P->col_coupleBlock;
711       Paso_Pattern *main_pattern=main_block->pattern;
712       Paso_Pattern *couple_pattern=couple_block->pattern;
713       const dim_t my_n=A->mainBlock->numRows;
714       index_t range, range_j;
715    
716     dim_t i, q;     dim_t i, q;
    const dim_t n =A_p->numRows;  
717     double *D_s=NULL;     double *D_s=NULL;
718     index_t *D_s_offset=NULL, iPtr, iPtr_j;     index_t *D_s_offset=NULL, iPtr, iPtr_j;
719     const dim_t ll = Paso_Util_iMax(n, degree_S);     const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A->mainBlock);
720     const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);     index_t *start_p_main_i, *start_p_couple_i;
721         dim_t degree_p_main_i, degree_p_couple_i;
722       dim_t len, ll, main_len, len_D_s;
723    
724       len = A->col_coupleBlock->numCols;
725       len = MAX(A->remote_coupleBlock->numCols, len);
726       ll = len + my_n;
727       main_len = main_pattern->len;
728    
729     #pragma omp parallel  private(D_s, D_s_offset, iPtr, q, iPtr_j)     #pragma omp parallel private(D_s,D_s_offset,i,start_p_main_i,start_p_couple_i,degree_p_main_i,degree_p_couple_i,range,iPtr,q,range_j,iPtr_j,len_D_s)
730     {     {
731          D_s=TMPMEMALLOC(ll,double);          D_s=TMPMEMALLOC(ll,double);
732          D_s_offset=TMPMEMALLOC(ll,index_t);          D_s_offset=TMPMEMALLOC(ll,index_t);
733    
734        #pragma omp for schedule(static)
735      #pragma omp for private(i) schedule(static)          for (i=0;i<my_n;++i) {
         for (i=0;i<n;++i) {  
736              if (counter_C[i]>=0) {              if (counter_C[i]>=0) {
737              P_p->val[P_p->pattern->ptr[i]]=1.;  /* i is a C row */              main_block->val[main_pattern->ptr[i]]=1.;  /* i is a C row */
738              } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {              } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
739             const index_t *start_s = &(S[offset_S[i]]);                 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
740             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);            /* this loop sums up the weak connections in a and creates
741                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];           a list of the strong connected columns which are not in
742                /* this loop sums up the weak connections in A and creates a list of the strong connected columns           C (=no interpolation nodes) */
743                   which are not in C (=no interpolation nodes) */            const index_t *start_s = &(S[offset_S[i]]);
744                const double A_ii = A_p->val[ptr_main_A[i]];                const double A_ii = A->mainBlock->val[ptr_main_A[i]];
745                double a=A_ii;                double a=A_ii;
746    
747            for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {            start_p_main_i = &(main_pattern->index[main_pattern->ptr[i]]);
748               const index_t j=A_p->pattern->index[iPtr];                start_p_couple_i = &(couple_pattern->index[couple_pattern->ptr[i]]);
749               const double A_ij=A_p->val[iPtr];                degree_p_main_i = main_pattern->ptr[i+1] - main_pattern->ptr[i];
750                  degree_p_couple_i = couple_pattern->ptr[i+1] - couple_pattern->ptr[i];
751    
752              /* first, check the mainBlock */
753              range = A->mainBlock->pattern->ptr[i + 1];
754              for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
755                 const index_t j=A->mainBlock->pattern->index[iPtr];
756                 const double A_ij=A->mainBlock->val[iPtr];
757                   if ( (i!=j) && (degree_S[j]>0) ) {                   if ( (i!=j) && (degree_S[j]>0) ) {
758                      /* is (i,j) a strong connection ?*/                      /* is (i,j) a strong connection ?*/
759                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);                  const index_t *where_s=(index_t*)bsearch(&j, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
# Line 408  void Paso_Preconditioner_AMG_setClassicP Line 761  void Paso_Preconditioner_AMG_setClassicP
761                          a+=A_ij;                            a+=A_ij;  
762                      } else {   /* yes i strongly connected with j */                      } else {   /* yes i strongly connected with j */
763                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */
764                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
765                                 if (where_p == NULL)  {                                 if (where_p == NULL)  {
766                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_setClassicProlongation: interpolation point is missing.");
767                                 } else {                                 } else {
768                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                              const index_t offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p_main_i);
769                                  P_p->val[offset]+=A_ij;                                  main_block->val[offset]+=A_ij;
770                                 }                                 }
771                            } else {  /* j is not an interpolation point */                            } else {  /* j is not an interpolation point */
772                                 /* find all interpolation points m of k */                                 /* find all interpolation points m of k */
                            const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]);  
                                const dim_t degree_P_j   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];  
773                                 double s=0.;                                 double s=0.;
774                                 dim_t len_D_s=0;                                 len_D_s=0;
775                             for (iPtr_j=A_p->pattern->ptr[j];iPtr_j<A_p->pattern->ptr[j + 1]; ++iPtr_j) {  
776                                  const double A_jm=A_p->val[iPtr_j];                     /* first, the mainBlock part */
777                                  const index_t m=A_p->pattern->index[iPtr_j];                     range_j = A->mainBlock->pattern->ptr[j + 1];
778                               for (iPtr_j=A->mainBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
779                                    const double A_jm=A->mainBlock->val[iPtr_j];
780                                    const index_t m=A->mainBlock->pattern->index[iPtr_j];
781                                      /* is m an interpolation point ? */                                      /* is m an interpolation point ? */
782                                  const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex);                                  const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
783                                      if (! (where_p_m==NULL)) {                                      if (! (where_p_m==NULL)) {
784                                   const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j);                                   const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
785                                           if (! SAMESIGN(A_ii,A_jm)) {                                           if (! SAMESIGN(A_ii,A_jm)) {
786                                                D_s[len_D_s]=A_jm;                                                D_s[len_D_s]=A_jm;
787                                           } else {                                           } else {
# Line 437  void Paso_Preconditioner_AMG_setClassicP Line 791  void Paso_Preconditioner_AMG_setClassicP
791                                           len_D_s++;                                           len_D_s++;
792                                      }                                      }
793                                 }                                 }
794    
795                       /* then the coupleBlock part */
796                       if (degree_p_couple_i) {
797                     range_j = A->col_coupleBlock->pattern->ptr[j + 1];
798                     for (iPtr_j=A->col_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
799                                        const double A_jm=A->col_coupleBlock->val[iPtr_j];
800                                        const index_t m=A->col_coupleBlock->pattern->index[iPtr_j];
801                                        /* is m an interpolation point ? */
802                                        const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
803                                        if (! (where_p_m==NULL)) {
804                                             const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
805                                             if (! SAMESIGN(A_ii,A_jm)) {
806                                                  D_s[len_D_s]=A_jm;
807                                             } else {
808                                                  D_s[len_D_s]=0.;
809                                             }
810                                             D_s_offset[len_D_s]=offset_m + main_len;
811                                             len_D_s++;
812                                        }
813                     }
814                                   }
815    
816                                 for (q=0;q<len_D_s;++q) s+=D_s[q];                                 for (q=0;q<len_D_s;++q) s+=D_s[q];
817                                 if (ABS(s)>0) {                                 if (ABS(s)>0) {
818                                     s=A_ij/s;                                     s=A_ij/s;
819                                     for (q=0;q<len_D_s;++q) {                                     for (q=0;q<len_D_s;++q) {
820                                          P_p->val[D_s_offset[q]]+=s*D_s[q];                      if (D_s_offset[q] < main_len)
821                                              main_block->val[D_s_offset[q]]+=s*D_s[q];
822                        else
823                          couple_block->val[D_s_offset[q]-main_len]+=s*D_s[q];
824                                     }                                     }
825                                 } else {                                 } else {
826                                     a+=A_ij;                                     a+=A_ij;
# Line 449  void Paso_Preconditioner_AMG_setClassicP Line 828  void Paso_Preconditioner_AMG_setClassicP
828                            }                            }
829                       }                       }
830                   }                   }
831                }  /* i has been processed, now we need to do some rescaling */                }  
832    
833              if (A->mpi_info->size > 1) {
834                  /* now, deal with the coupleBlock */
835                  range = A->col_coupleBlock->pattern->ptr[i + 1];
836                  for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
837                     const index_t j=A->col_coupleBlock->pattern->index[iPtr];
838                     const double A_ij=A->col_coupleBlock->val[iPtr];
839                     if ( (i!=j) && (degree_S[j]>0) ) {
840                        /* is (i,j) a strong connection ?*/
841                index_t t=j+my_n;
842                        const index_t *where_s=(index_t*)bsearch(&t, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
843                        if (where_s == NULL) { /* weak connections are accummulated */
844                            a+=A_ij;
845                        } else {   /* yes i stronly connect with j */
846                            if  (counter_C[t]>=0)  { /* j is an interpolation point : add A_ij into P */
847                                   const index_t *where_p=(index_t*)bsearch(&counter_C[t], start_p_couple_i,degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
848                                   if (where_p == NULL)  {
849                                           Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");
850                                   } else {
851                                        const index_t offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p_couple_i);
852                                        couple_block->val[offset]+=A_ij;
853                                   }
854                              } else {  /* j is not an interpolation point */
855                                   /* find all interpolation points m of k */
856                       double s=0.;
857                                   len_D_s=0;
858    
859                                   /* first, the row_coupleBlock part */
860                                   range_j = A->row_coupleBlock->pattern->ptr[j + 1];
861                                   for (iPtr_j=A->row_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
862                                        const double A_jm=A->row_coupleBlock->val[iPtr_j];
863                                        const index_t m=A->row_coupleBlock->pattern->index[iPtr_j];
864                                        /* is m an interpolation point ? */
865                                        const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
866                                        if (! (where_p_m==NULL)) {
867                                             const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
868                                             if (! SAMESIGN(A_ii,A_jm)) {
869                                                  D_s[len_D_s]=A_jm;
870                                             } else {
871                                                  D_s[len_D_s]=0.;
872                                             }
873                                             D_s_offset[len_D_s]=offset_m;
874                                             len_D_s++;
875                                        }
876                                   }
877    
878                                   /* then the remote_coupleBlock part */
879                       range_j = A->remote_coupleBlock->pattern->ptr[j + 1];
880                       for (iPtr_j=A->remote_coupleBlock->pattern->ptr[j]; iPtr_j<range_j; iPtr_j++) {
881                                        const double A_jm=A->remote_coupleBlock->val[iPtr_j];
882                                        const index_t m=A->remote_coupleBlock->pattern->index[iPtr_j];
883                                        /* is m an interpolation point ? */
884                                        const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i, degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
885                                        if (! (where_p_m==NULL)) {
886                                             const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
887                                             if (! SAMESIGN(A_ii,A_jm)) {
888                                                  D_s[len_D_s]=A_jm;
889                                             } else {
890                                                  D_s[len_D_s]=0.;
891                                             }
892                                             D_s_offset[len_D_s]=offset_m + main_len;
893                                             len_D_s++;
894                                        }
895                                   }
896    
897                                   for (q=0;q<len_D_s;++q) s+=D_s[q];
898                                   if (ABS(s)>0) {
899                                       s=A_ij/s;
900                                       for (q=0;q<len_D_s;++q) {
901                        if (D_s_offset[q] < main_len)
902                                              main_block->val[D_s_offset[q]]+=s*D_s[q];
903                        else
904                          couple_block->val[D_s_offset[q]-main_len]+=s*D_s[q];
905                                       }
906                                   } else {
907                                       a+=A_ij;
908                                   }
909                              }
910                         }
911                     }
912                  }
913              }
914    
915              /* i has been processed, now we need to do some rescaling */
916                if (ABS(a)>0.) {                if (ABS(a)>0.) {
917                     a=-1./a;                     a=-1./a;
918                     for (iPtr=P_p->pattern->ptr[i]; iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {             range = main_pattern->ptr[i + 1];
919                          P_p->val[iPtr]*=a;                     for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
920                            main_block->val[iPtr]*=a;
921                       }
922    
923               range = couple_pattern->ptr[i + 1];
924                       for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
925                            couple_block->val[iPtr]*=a;
926                     }                     }
927                }                }
928            }            }
# Line 463  void Paso_Preconditioner_AMG_setClassicP Line 932  void Paso_Preconditioner_AMG_setClassicP
932       }    /* end of parallel region */       }    /* end of parallel region */
933  }  }
934    
935  void Paso_Preconditioner_AMG_setClassicProlongation_Block(Paso_SparseMatrix* P_p,  void Paso_Preconditioner_AMG_setClassicProlongation_Block(
936                                 Paso_SparseMatrix* A_p,      Paso_SystemMatrix* P, Paso_SystemMatrix* A, const index_t* offset_S,
937                                                     const index_t* offset_S, const dim_t* degree_S, const index_t* S,      const dim_t* degree_S, const index_t* S, const index_t *counter_C) {
938                             const index_t *counter_C) {     Paso_SparseMatrix *main_block=P->mainBlock;
939       Paso_SparseMatrix *couple_block=P->col_coupleBlock;
940       Paso_Pattern *main_pattern=main_block->pattern;
941       Paso_Pattern *couple_pattern=couple_block->pattern;
942       const dim_t row_block=A->row_block_size;
943       const dim_t my_n=A->mainBlock->numRows;
944       const dim_t A_block = A->block_size;
945       index_t range, range_j;
946    
947     dim_t i, q, ib;     dim_t i, q, ib;
    const dim_t row_block=A_p->row_block_size;  
    const dim_t A_block = A_p->block_size;  
    const dim_t n =A_p->numRows;  
948     double *D_s=NULL;     double *D_s=NULL;
949     index_t *D_s_offset=NULL, iPtr, iPtr_j;     index_t *D_s_offset=NULL, iPtr, iPtr_j;
950     const dim_t ll = Paso_Util_iMax(n, degree_S);     index_t *start_p_main_i, *start_p_couple_i;
951     const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A_p);     dim_t degree_p_main_i, degree_p_couple_i;
952       dim_t len, ll, main_len, len_D_s;
953       const index_t *ptr_main_A = Paso_SparseMatrix_borrowMainDiagonalPointer(A->mainBlock);
954        
955       len = A->col_coupleBlock->numCols;
956     #pragma omp parallel  private(D_s, D_s_offset, iPtr, q, iPtr_j,ib)     len = MAX(A->remote_coupleBlock->numCols, len);
957       ll = len + my_n;
958       main_len = main_pattern->len;
959       #pragma omp parallel private(D_s,D_s_offset,i,ib,start_p_main_i,start_p_couple_i,degree_p_main_i,degree_p_couple_i,range,iPtr,q,range_j,iPtr_j,len_D_s)
960     {     {
961          double *a=TMPMEMALLOC(row_block, double);          double *a=TMPMEMALLOC(row_block, double);
962          D_s=TMPMEMALLOC(row_block*ll,double);          D_s=TMPMEMALLOC(row_block*ll,double);
963          D_s_offset=TMPMEMALLOC(row_block*ll,index_t);          D_s_offset=TMPMEMALLOC(row_block*ll,index_t);
964    
965      #pragma omp for private(i) schedule(static)      #pragma omp for private(i) schedule(static)
966          for (i=0;i<n;++i) {          for (i=0;i<my_n;++i) {
967              if (counter_C[i]>=0) {              if (counter_C[i]>=0) {
968              const index_t offset = P_p->pattern->ptr[i];              const index_t offset = main_pattern->ptr[i];
969              for (ib =0; ib<row_block; ++ib) P_p->val[row_block*offset+ib]=1.;  /* i is a C row */              for (ib =0; ib<row_block; ++ib) main_block->val[row_block*offset+ib]=1.;  /* i is a C row */
970              } else if (P_p->pattern->ptr[i + 1] > P_p->pattern->ptr[i]) {              } else if ((main_pattern->ptr[i + 1] > main_pattern->ptr[i]) ||
971             const index_t *start_s = &(S[offset_S[i]]);                 (couple_pattern->ptr[i + 1] > couple_pattern->ptr[i])) {
972             const index_t *start_p = &(P_p->pattern->index[P_p->pattern->ptr[i]]);            /* this loop sums up the weak connections in a and creates
973                 const dim_t degree_P_i   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];           a list of the strong connected columns which are not in
974                /* this loop sums up the weak connections in A and creates a list of the strong connected columns           C (=no interpolation nodes) */
975                   which are not in C (=no interpolation nodes) */            const index_t *start_s = &(S[offset_S[i]]);
976                const double *A_ii = &(A_p->val[ptr_main_A[i]*A_block]);            const double *A_ii = &(A->mainBlock->val[ptr_main_A[i]*A_block]);
977              start_p_main_i = &(main_pattern->index[main_pattern->ptr[i]]);
978              start_p_couple_i = &(couple_pattern->index[couple_pattern->ptr[i]]);
979                  degree_p_main_i = main_pattern->ptr[i+1] - main_pattern->ptr[i];
980              degree_p_couple_i = couple_pattern->ptr[i+1] - couple_pattern->ptr[i];
981                for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];                for (ib=0; ib<row_block; ib++) a[ib]=A_ii[(row_block+1)*ib];
                 
982    
983            for (iPtr=A_p->pattern->ptr[i];iPtr<A_p->pattern->ptr[i + 1]; ++iPtr) {            /* first, check the mainBlock */
984               const index_t j=A_p->pattern->index[iPtr];            range = A->mainBlock->pattern->ptr[i + 1];
985               const double* A_ij=&(A_p->val[iPtr*A_block]);            for (iPtr=A->mainBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
986                 const index_t j=A->mainBlock->pattern->index[iPtr];
987                 const double* A_ij=&(A->mainBlock->val[iPtr*A_block]);
988    
989                   if ( (i!=j) && (degree_S[j]>0) ) {                   if ( (i!=j) && (degree_S[j]>0) ) {
990                      /* is (i,j) a strong connection ?*/                      /* is (i,j) a strong connection ?*/
# Line 509  void Paso_Preconditioner_AMG_setClassicP Line 993  void Paso_Preconditioner_AMG_setClassicP
993                          for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];                          for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
994                      } else {   /* yes i strongly connected with j */                      } else {   /* yes i strongly connected with j */
995                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */                          if  (counter_C[j]>=0)  { /* j is an interpolation point : add A_ij into P */
996                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p,degree_P_i, sizeof(index_t), Paso_comparIndex);                             const index_t *where_p=(index_t*)bsearch(&counter_C[j], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
997                                 if (where_p == NULL)  {                                 if (where_p == NULL)  {
998                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation: interpolation point is missing.");                                         Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
999                                 } else {                                 } else {
1000                              const index_t offset = P_p->pattern->ptr[i]+ (index_t)(where_p-start_p);                              const index_t offset = main_pattern->ptr[i]+ (index_t)(where_p-start_p_main_i);
1001                                      for (ib=0; ib<row_block; ib++) P_p->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];                                      for (ib=0; ib<row_block; ib++) main_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
1002                                 }                                 }
1003                            } else {  /* j is not an interpolation point */                            } else {  /* j is not an interpolation point */
1004                                 /* find all interpolation points m of k */                                 /* find all interpolation points m of k */
1005                             const index_t *start_p_j = &(P_p->pattern->index[P_p->pattern->ptr[i]]);                                 len_D_s=0;
1006                                 const dim_t degree_P_j   = P_p->pattern->ptr[i + 1]-P_p->pattern->ptr[i];  
1007                                 dim_t len_D_s=0;                     /* first, the mainBlock part */
1008                             for (iPtr_j=A_p->pattern->ptr[j];iPtr_j<A_p->pattern->ptr[j + 1]; ++iPtr_j) {                     range_j = A->mainBlock->pattern->ptr[j + 1];
1009                                  const double* A_jm=&(A_p->val[iPtr_j*A_block]);                             for (iPtr_j=A->mainBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1010                                  const index_t m=A_p->pattern->index[iPtr_j];                                  const double* A_jm=&(A->mainBlock->val[iPtr_j*A_block]);
1011                                    const index_t m=A->mainBlock->pattern->index[iPtr_j];
1012                                        /* is m an interpolation point ? */
1013                                    const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
1014                                        if (! (where_p_m==NULL)) {
1015                                     const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1016                                             for (ib=0; ib<row_block; ib++) {
1017                                                  if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1018                                                       D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1019                                                  } else {
1020                                                       D_s[len_D_s*row_block+ib]=0.;
1021                                                  }
1022                                             }
1023                                             D_s_offset[len_D_s] = offset_m;
1024                                             len_D_s++;
1025                                        }
1026                                   }
1027    
1028                       /* then the coupleBlock part */
1029                       range_j = A->col_coupleBlock->pattern->ptr[j+1];
1030                       for (iPtr_j=A->col_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1031                                        const double* A_jm=&(A->col_coupleBlock->val[iPtr_j*A_block]);
1032                                        const index_t m=A->col_coupleBlock->pattern->index[iPtr_j];
1033                                      /* is m an interpolation point ? */                                      /* is m an interpolation point ? */
1034                                  const index_t *where_p_m=(index_t*)bsearch(&counter_C[m], start_p_j,degree_P_j, sizeof(index_t), Paso_comparIndex);                                      const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
1035                                      if (! (where_p_m==NULL)) {                                      if (! (where_p_m==NULL)) {
1036                                   const index_t offset_m = P_p->pattern->ptr[i]+ (index_t)(where_p_m-start_p_j);                                           const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1037                                             for (ib=0; ib<row_block; ib++) {
1038                                                  if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1039                                                       D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1040                                                  } else {
1041                                                       D_s[len_D_s*row_block+ib]=0.;
1042                                                  }
1043                                             }
1044                                             D_s_offset[len_D_s]=offset_m + main_len;
1045                                             len_D_s++;
1046                                        }
1047                       }
1048    
1049                                   for (ib=0; ib<row_block; ib++) {
1050                                       double s=0;
1051                                       for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1052                            
1053                                       if (ABS(s)>0) {
1054                                         s=A_ij[(row_block+1)*ib]/s;
1055                         for (q=0; q<len_D_s; q++) {
1056                           if (D_s_offset[q] < main_len)
1057                                                main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1058                           else{
1059                            couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1060                        }
1061                         }
1062                                       } else {
1063                                           a[ib]+=A_ij[(row_block+1)*ib];
1064                                       }
1065                                   }
1066                              }
1067                         }
1068                     }
1069                  }  
1070    
1071              if (A->mpi_info->size > 1) {
1072              /* now, deal with the coupleBlock */
1073              range = A->col_coupleBlock->pattern->ptr[i + 1];
1074              for (iPtr=A->col_coupleBlock->pattern->ptr[i]; iPtr<range; iPtr++) {
1075                 const index_t j=A->col_coupleBlock->pattern->index[iPtr];
1076                 const double* A_ij=&(A->col_coupleBlock->val[iPtr*A_block]);
1077    
1078                     if ( (i!=j) && (degree_S[j]>0) ) {
1079                        /* is (i,j) a strong connection ?*/
1080                index_t t=j+my_n;
1081                    const index_t *where_s=(index_t*)bsearch(&t, start_s,degree_S[i],sizeof(index_t), Paso_comparIndex);
1082                    if (where_s == NULL) { /* weak connections are accummulated */
1083                            for (ib=0; ib<row_block; ib++) a[ib]+=A_ij[(row_block+1)*ib];
1084                        } else {   /* yes i strongly connected with j */
1085                            if  (counter_C[t]>=0)  { /* j is an interpolation point : add A_ij into P */
1086                               const index_t *where_p=(index_t*)bsearch(&counter_C[t], start_p_couple_i,degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
1087                                   if (where_p == NULL)  {
1088                                           Esys_setError(SYSTEM_ERROR, "Paso_Preconditioner_AMG_setClassicProlongation_Block: interpolation point is missing.");
1089    
1090                                   } else {
1091                                const index_t offset = couple_pattern->ptr[i]+ (index_t)(where_p-start_p_couple_i);
1092                                        for (ib=0; ib<row_block; ib++) couple_block->val[offset*row_block+ib] +=A_ij[(row_block+1)*ib];
1093                                   }
1094                              } else {  /* j is not an interpolation point */
1095                                   /* find all interpolation points m of k */
1096                                   len_D_s=0;
1097    
1098                       /* first, the row_coupleBlock part */
1099                       range_j = A->row_coupleBlock->pattern->ptr[j + 1];
1100                               for (iPtr_j=A->row_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1101                                        /* is m an interpolation point ? */
1102                        index_t m, *where_p_m;
1103                        double *A_jm;
1104                        A_jm=&(A->row_coupleBlock->val[iPtr_j*A_block]);
1105                        m=A->row_coupleBlock->pattern->index[iPtr_j];
1106    
1107                        where_p_m=(index_t*)bsearch(&counter_C[m], start_p_main_i,degree_p_main_i, sizeof(index_t), Paso_comparIndex);
1108                                        if (! (where_p_m==NULL)) {
1109                                     const index_t offset_m = main_pattern->ptr[i]+ (index_t)(where_p_m-start_p_main_i);
1110                                           for (ib=0; ib<row_block; ib++) {                                           for (ib=0; ib<row_block; ib++) {
1111                                                if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {                                                if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1112                                                     D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];                                                     D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
# Line 539  void Paso_Preconditioner_AMG_setClassicP Line 1118  void Paso_Preconditioner_AMG_setClassicP
1118                                           len_D_s++;                                           len_D_s++;
1119                                      }                                      }
1120                                 }                                 }
1121    
1122                       /* then the remote_coupleBlock part */
1123                       if (degree_p_couple_i) {
1124                     range_j = A->remote_coupleBlock->pattern->ptr[j + 1];
1125                     for (iPtr_j=A->remote_coupleBlock->pattern->ptr[j];iPtr_j<range_j; iPtr_j++) {
1126                                        const double* A_jm=&(A->remote_coupleBlock->val[iPtr_j*A_block]);
1127                                        const index_t m=A->remote_coupleBlock->pattern->index[iPtr_j];
1128                                        /* is m an interpolation point ? */
1129                                        const index_t *where_p_m=(index_t*)bsearch(&counter_C[m+my_n], start_p_couple_i,degree_p_couple_i, sizeof(index_t), Paso_comparIndex);
1130                                        if (! (where_p_m==NULL)) {
1131                                             const index_t offset_m = couple_pattern->ptr[i]+ (index_t)(where_p_m-start_p_couple_i);
1132                                             for (ib=0; ib<row_block; ib++) {
1133                                                  if (! SAMESIGN(A_ii[(row_block+1)*ib],A_jm[(row_block+1)*ib]) ) {
1134                                                       D_s[len_D_s*row_block+ib]=A_jm[(row_block+1)*ib];
1135                                                  } else {
1136                                                       D_s[len_D_s*row_block+ib]=0.;
1137                                                  }
1138                                             }
1139                                             D_s_offset[len_D_s]=offset_m + main_len;
1140                                             len_D_s++;
1141                                        }
1142                     }
1143                                   }
1144    
1145                                 for (ib=0; ib<row_block; ib++) {                                 for (ib=0; ib<row_block; ib++) {
1146                                     double s=0;                                     double s=0;
1147                                     for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];                                     for (q=0;q<len_D_s;++q) s+=D_s[q*row_block+ib];
1148                                                    
1149                                     if (ABS(s)>0) {                                     if (ABS(s)>0) {
1150                                         s=A_ij[(row_block+1)*ib]/s;                                         s=A_ij[(row_block+1)*ib]/s;
1151                                         for (q=0;q<len_D_s;++q) {                         for (q=0;q<len_D_s;++q) {
1152                                              P_p->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];                       if (D_s_offset[q] < main_len)
1153                                         }                                              main_block->val[D_s_offset[q]*row_block+ib]+=s*D_s[q*row_block+ib];
1154                         else
1155                            couple_block->val[(D_s_offset[q]-main_len)*row_block+ib]+=s*D_s[q*row_block+ib];
1156                           }
1157                                     } else {                                     } else {
1158                                         a[ib]+=A_ij[(row_block+1)*ib];                                         a[ib]+=A_ij[(row_block+1)*ib];
1159                                     }                                     }
# Line 555  void Paso_Preconditioner_AMG_setClassicP Line 1161  void Paso_Preconditioner_AMG_setClassicP
1161                            }                            }
1162                       }                       }
1163                   }                   }
1164                }  /* i has been processed, now we need to do some rescaling */            }
1165              }
1166    
1167              /* i has been processed, now we need to do some rescaling */
1168                for (ib=0; ib<row_block; ib++) {                for (ib=0; ib<row_block; ib++) {
1169                     register double a2=a[ib];                     register double a2=a[ib];
1170                     if (ABS(a2)>0.) {                     if (ABS(a2)>0.) {
1171                          a2=-1./a2;                          a2=-1./a2;
1172                          for (iPtr=P_p->pattern->ptr[i]; iPtr<P_p->pattern->ptr[i + 1]; ++iPtr) {              range = main_pattern->ptr[i + 1];
1173                               P_p->val[iPtr*row_block+ib]*=a2;                          for (iPtr=main_pattern->ptr[i]; iPtr<range; iPtr++) {
1174                                 main_block->val[iPtr*row_block+ib]*=a2;
1175                            }
1176    
1177                range = couple_pattern->ptr[i + 1];
1178                for (iPtr=couple_pattern->ptr[i]; iPtr<range; iPtr++) {
1179                                 couple_block->val[iPtr*row_block+ib]*=a2;
1180                          }                          }
1181                     }                     }
1182                }                }

Legend:
Removed from v.3867  
changed lines
  Added in v.3868

  ViewVC Help
Powered by ViewVC 1.1.26