/[escript]/trunk/paso/src/Smoother.c
ViewVC logotype

Diff of /trunk/paso/src/Smoother.c

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

revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2998 by artak, Wed Mar 24 04:00:20 2010 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 77  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 77  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
77    
78            if (! (Paso_checkPtr(out->diag))) {            if (! (Paso_checkPtr(out->diag))) {
79               if (n_block==1) {               if (n_block==1) {
80                  #pragma omp parallel for private(i, iPtr) schedule(static)                  #pragma omp parallel for private(i,iPtr,iptr_main) schedule(static)
81                  for (i = 0; i < A->pattern->numOutput; i++) {                  for (i = 0; i < A->pattern->numOutput; i++) {
82                       iptr_main=0;
83                     out->diag[i]=1.;                     out->diag[i]=1.;
84                     /* find main diagonal */                     /* find main diagonal */
85                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {
# Line 91  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 92  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
92                     out->main_iptr[i]=iptr_main;                     out->main_iptr[i]=iptr_main;
93                  }                  }
94               } else if (n_block==2) {               } else if (n_block==2) {
95                  #pragma omp parallel for private(i, iPtr) schedule(static)                  #pragma omp parallel for private(i,iPtr,iptr_main) schedule(static)
96                  for (i = 0; i < A->pattern->numOutput; i++) {                  for (i = 0; i < A->pattern->numOutput; i++) {
97                     out->diag[i*4+0]= 1.;                     out->diag[i*4+0]= 1.;
98                     out->diag[i*4+1]= 0.;                     out->diag[i*4+1]= 0.;
99                     out->diag[i*4+2]= 0.;                     out->diag[i*4+2]= 0.;
100                     out->diag[i*4+3]= 1.;                     out->diag[i*4+3]= 1.;
101                       iptr_main=0;
102                     /* find main diagonal */                     /* find main diagonal */
103                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {
104                         if (A->pattern->index[iPtr]==i) {                         if (A->pattern->index[iPtr]==i) {
# Line 111  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 113  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
113                     out->main_iptr[i]=iptr_main;                     out->main_iptr[i]=iptr_main;
114                  }                    }  
115               } else if (n_block==3) {               } else if (n_block==3) {
116                  #pragma omp parallel for private(i, iPtr) schedule(static)                  #pragma omp parallel for private(i, iPtr,iptr_main) schedule(static)
117                  for (i = 0; i < A->pattern->numOutput; i++) {                  for (i = 0; i < A->pattern->numOutput; i++) {
118                     out->diag[i*9  ]=1.;                     out->diag[i*9  ]=1.;
119                     out->diag[i*9+1]=0.;                     out->diag[i*9+1]=0.;
# Line 122  Paso_Solver_GS* Paso_Solver_getGS(Paso_S Line 124  Paso_Solver_GS* Paso_Solver_getGS(Paso_S
124                     out->diag[i*9+6]=0.;                     out->diag[i*9+6]=0.;
125                     out->diag[i*9+7]=0.;                     out->diag[i*9+7]=0.;
126                     out->diag[i*9+8]=1.;                     out->diag[i*9+8]=1.;
127                       iptr_main=0;
128                     /* find main diagonal */                     /* find main diagonal */
129                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {                     for (iPtr = A->pattern->ptr[i]; iPtr < A->pattern->ptr[i + 1]; iPtr++) {
130                         if (A->pattern->index[iPtr]==i) {                         if (A->pattern->index[iPtr]==i) {
# Line 200  void Paso_Solver_solveGS(Paso_Solver_GS Line 203  void Paso_Solver_solveGS(Paso_Solver_GS
203                     }                     }
204                }                }
205             } else if (n_block==2) {             } else if (n_block==2) {
206                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,iptr_main,S1,S2,R1,R2,A11,A21,A12,A22,D,S11,S21,S12,S22)
207                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
208                     if (gs->colorOf[i]==color) {                     if (gs->colorOf[i]==color) {
209                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 236  void Paso_Solver_solveGS(Paso_Solver_GS Line 239  void Paso_Solver_solveGS(Paso_Solver_GS
239    
240                }                }
241             } else if (n_block==3) {             } else if (n_block==3) {
242                #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3)                #pragma omp parallel for schedule(static) private(i,iptr_ik,iptr_main,k,S1,S2,S3,R1,R2,R3,A11,A21,A31,A12,A22,A32,A13,A23,A33,D,S11,S21,S31,S12,S22,S32,S13,S23,S33)
243                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
244                     if (gs->colorOf[i]==color) {                     if (gs->colorOf[i]==color) {
245                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 294  void Paso_Solver_solveGS(Paso_Solver_GS Line 297  void Paso_Solver_solveGS(Paso_Solver_GS
297       /* backward substitution */       /* backward substitution */
298       for (color=(gs->num_colors)-1;color>-1;--color) {       for (color=(gs->num_colors)-1;color>-1;--color) {
299             if (n_block==1) {             if (n_block==1) {
300                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1)                /*#pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,R1,iptr_main)*/
301                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
302                     if (gs->colorOf[i]==color) {                     if (gs->colorOf[i]==color) {
303                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 312  void Paso_Solver_solveGS(Paso_Solver_GS Line 315  void Paso_Solver_solveGS(Paso_Solver_GS
315                     }                     }
316                }                }
317             } else if (n_block==2) {             } else if (n_block==2) {
318                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,R1,R2,iptr_main,D,A11,A21,A12,A22,S11,S21,S12,S22)
319                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
320                     if (gs->colorOf[i]==color) {                     if (gs->colorOf[i]==color) {
321                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */
# Line 350  void Paso_Solver_solveGS(Paso_Solver_GS Line 353  void Paso_Solver_solveGS(Paso_Solver_GS
353                      }                      }
354                }                }
355             } else if (n_block==3) {             } else if (n_block==3) {
356                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3)                #pragma omp parallel for schedule(static) private(i,iptr_ik,k,S1,S2,S3,R1,R2,R3,iptr_main,D,A11,A21,A31,A12,A22,A32,A13,A23,A33,S11,S21,S31,S12,S22,S32,S13,S23,S33)
357                for (i = 0; i < n; ++i) {                for (i = 0; i < n; ++i) {
358                     if (gs->colorOf[i]==color) {                     if (gs->colorOf[i]==color) {
359                       /* x_i=x_i-a_ik*x_k */                       /* x_i=x_i-a_ik*x_k */

Legend:
Removed from v.2548  
changed lines
  Added in v.2998

  ViewVC Help
Powered by ViewVC 1.1.26