/[escript]/branches/stage3.1/paso/src/SolverFCT_FluxControl.c
ViewVC logotype

Diff of /branches/stage3.1/paso/src/SolverFCT_FluxControl.c

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

revision 2944 by jfenwick, Thu Feb 4 01:42:47 2010 UTC revision 2945 by jfenwick, Wed Feb 24 00:17:46 2010 UTC
# Line 64  void Paso_FCTransportProblem_setLowOrder Line 64  void Paso_FCTransportProblem_setLowOrder
64        for (i = 0; i < n; ++i) {        for (i = 0; i < n; ++i) {
65            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];            sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
66            /* look at a[i,j] */            /* look at a[i,j] */
           #pragma ivdep  
67            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {            for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
68                j=pattern->mainPattern->index[iptr_ij];                j=pattern->mainPattern->index[iptr_ij];
69                if (j!=i) {                if (j!=i) {
70                   /* find entry a[j,i] */                   /* find entry a[j,i] */
71                     #pragma ivdep
72                   for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {                   for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
73                      if (pattern->mainPattern->index[iptr_ji]==i) {                      if ( pattern->mainPattern->index[iptr_ji] == i) {
74                          rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];                          rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
75                          rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];                          rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
76                          d_ij=-MIN3(0.,rtmp1,rtmp2);                          d_ij=-MIN3(0.,rtmp1,rtmp2);
# Line 80  void Paso_FCTransportProblem_setLowOrder Line 80  void Paso_FCTransportProblem_setLowOrder
80                      }                      }
81                   }                   }
82               }               }
83           }            }
           #pragma ivdep  
84            for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {            for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
85                j=pattern->col_couplePattern->index[iptr_ij];                j=pattern->col_couplePattern->index[iptr_ij];
86                /* find entry a[j,i] */                /* find entry a[j,i] */
87                  #pragma ivdep
88                for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {                for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
89                      if (pattern->row_couplePattern->index[iptr_ji]==i) {                      if (pattern->row_couplePattern->index[iptr_ji]==i) {
90                          rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];                          rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
# Line 103  void Paso_FCTransportProblem_setLowOrder Line 103  void Paso_FCTransportProblem_setLowOrder
103    }    }
104  }  }
105  /*  /*
106   * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i   * out_i=m_i u_i + a * \sum_{j <> i} l_{ij} (u_j-u_i) + b * q_i
107   *   *
108   */   */
109  void Paso_SolverFCT_setMuPaLuPbQ(double* out,  void Paso_SolverFCT_setMuPaLuPbQ(double* out,
# Line 120  void Paso_SolverFCT_setMuPaLuPbQ(double* Line 120  void Paso_SolverFCT_setMuPaLuPbQ(double*
120    const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);    const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
121    register double sum, u_i, l_ij;    register double sum, u_i, l_ij;
122    register index_t iptr_ij;    register index_t iptr_ij;
   
123    n=Paso_SystemMatrix_getTotalNumRows(L);    n=Paso_SystemMatrix_getTotalNumRows(L);
124    
125    #pragma omp parallel for private(i) schedule(static)    #pragma omp parallel for private(i) schedule(static)
# Line 409  void Paso_FCTransportProblem_addCorrecte Line 408  void Paso_FCTransportProblem_addCorrecte
408    dim_t n, i, j;    dim_t n, i, j;
409    Paso_SystemMatrixPattern *pattern;    Paso_SystemMatrixPattern *pattern;
410    register double RN_i, RP_i, f_i, f_ij;    register double RN_i, RP_i, f_i, f_ij;
411      double alpha_i;
412    register index_t iptr_ij;    register index_t iptr_ij;
413    const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);    const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
414    const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);    const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
# Line 416  void Paso_FCTransportProblem_addCorrecte Line 416  void Paso_FCTransportProblem_addCorrecte
416    const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);    const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
417    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);    n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
418    
419    pattern=flux_matrix->pattern;    pattern=flux_matrix->pattern;
420    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)    #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
421    for (i = 0; i < n; ++i) {    for (i = 0; i < n; ++i) {
422        
423         alpha_i=0;
424       RN_i=RN[i];       RN_i=RN[i];
425       RP_i=RP[i];       RP_i=RP[i];
426       f_i=0;       f_i=0;
# Line 428  void Paso_FCTransportProblem_addCorrecte Line 430  void Paso_FCTransportProblem_addCorrecte
430           f_ij=flux_matrix->mainBlock->val[iptr_ij];           f_ij=flux_matrix->mainBlock->val[iptr_ij];
431           if (f_ij >=0) {           if (f_ij >=0) {
432                f_i+=f_ij*MIN(RP_i,RN[j]);                f_i+=f_ij*MIN(RP_i,RN[j]);
433    /* printf("alpha %d %d = %e %e\n",i,j,MIN(RP_i,RN[j]),f_ij);          
434    alpha_i=MAX(alpha_i,MIN(RP_i,RN[j]) ); */
435           } else {           } else {
436                f_i+=f_ij*MIN(RN_i,RP[j]);                f_i+=f_ij*MIN(RN_i,RP[j]);
437           }  /*
438    printf("alpha %d %d = %e %e\n",i,j,MIN(RN_i,RP[j]),f_ij);  
439    alpha_i=MAX(alpha_i,MIN(RN_i,RP[j]) );
440    */
441         }
442       }       }
443       #pragma ivdep       #pragma ivdep
444       for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {       for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
# Line 442  void Paso_FCTransportProblem_addCorrecte Line 450  void Paso_FCTransportProblem_addCorrecte
450                f_i+=f_ij*MIN(RN_i,remote_RP[j]);                f_i+=f_ij*MIN(RN_i,remote_RP[j]);
451            }            }
452        }        }
453    /* printf("alpha %d = %e -> %e (%e, %e)\n",i,alpha_i, f_i, f[i],f[i]+f_i); */
454        f[i]+=f_i;        f[i]+=f_i;
455    }    }
456  }  }

Legend:
Removed from v.2944  
changed lines
  Added in v.2945

  ViewVC Help
Powered by ViewVC 1.1.26