/[escript]/trunk-mpi-branch/paso/src/SystemMatrix_MatrixVector.c
ViewVC logotype

Diff of /trunk-mpi-branch/paso/src/SystemMatrix_MatrixVector.c

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

revision 1164 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1165 by ksteube, Thu May 24 00:45:48 2007 UTC
# Line 497  void  Paso_SystemMatrix_MatrixVector_CSR Line 497  void  Paso_SystemMatrix_MatrixVector_CSR
497             index_t min_index, index_t max_index,             index_t min_index, index_t max_index,
498             double* out)             double* out)
499  {  {
500      register index_t ir,icol,iptr,icb,irb,irow,ic,Aiptr, itmp;      register index_t iptr,icb,irb,irow,ic,Aiptr,itmp,irtmp,ictmp;
501      register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;      register double reg,reg1,reg2,reg3,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22;
502      if (A ->col_block_size==1 && A->row_block_size ==1) {      if (A ->col_block_size==1 && A->row_block_size ==1) {
503          #pragma omp for private(irow,iptr,reg) schedule(static)      int iptr0;
504            #pragma omp for private(irow,iptr,reg,itmp,iptr0) schedule(static)
505      for (irow=0;irow< A->pattern->myNumOutput;++irow) {      for (irow=0;irow< A->pattern->myNumOutput;++irow) {
506            reg=0.;            reg=0.;
507            #pragma ivdep            #pragma ivdep
508        for (iptr=(A->pattern->ptr[irow]);iptr<(A->pattern->ptr[irow+1]); ++iptr) {        /* index[] is an increasing sequence, can find min first */
509          for (iptr0=(A->pattern->ptr[irow]);iptr0<(A->pattern->ptr[irow+1]); ++iptr0) {
510                  itmp=A->pattern->index[iptr0];
511                  if ( min_index<= itmp) { break; }
512          }
513              #pragma ivdep
514          for (iptr=iptr0;iptr<(A->pattern->ptr[irow+1]); ++iptr) {
515                itmp=A->pattern->index[iptr];                itmp=A->pattern->index[iptr];
516                if (( min_index<= itmp) && (itmp<max_index) ) {            /* index[] is an increasing sequence, can quit loop if at max */
517               reg += A->val[iptr] * in[itmp-min_index];            if (itmp >= max_index) { break; }
518                }            reg += A->val[iptr] * in[itmp-min_index];
519        }        }
520        out[irow] += alpha * reg;        out[irow] += alpha * reg;
521      }      }
522      } else if (A ->col_block_size==2 && A->row_block_size ==2) {      } else if (A ->col_block_size==2 && A->row_block_size ==2) {
523          #pragma omp for private(ir,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11) schedule(static)      int iptr0;
524      for (ir=0;ir< A->pattern->myNumOutput;ir++) {          #pragma omp for private(irow,reg1,reg2,iptr,ic,Aiptr,in1,in2,A00,A10,A01,A11,itmp,iptr0) schedule(static)
525        for (irow=0;irow< A->pattern->myNumOutput;++irow) {
526            reg1=0.;            reg1=0.;
527            reg2=0.;            reg2=0.;
528            #pragma ivdep            #pragma ivdep
529        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {        /* index[] is an increasing sequence, can find min first */
530          for (iptr0=(A->pattern->ptr[irow]);iptr0<(A->pattern->ptr[irow+1]); ++iptr0) {
531                  itmp=A->pattern->index[iptr0];
532                  if ( min_index<= itmp) { break; }
533          }
534              #pragma ivdep
535          for (iptr=iptr0;iptr<(A->pattern->ptr[irow+1]); ++iptr) {
536                itmp=A->pattern->index[iptr];                itmp=A->pattern->index[iptr];
537                if (( min_index<= itmp) && (itmp<max_index) ) {            /* index[] is an increasing sequence, can quit loop if at max */
538                ic=2*(itmp-min_index);            if (itmp >= max_index) { break; }
539                    Aiptr=iptr*4;            ic=2*(itmp-min_index);
540                    in1=in[ic];                Aiptr=iptr*4;
541                    in2=in[1+ic];                in1=in[ic];
542                    A00=A->val[Aiptr  ];                in2=in[1+ic];
543                    A10=A->val[Aiptr+1];                A00=A->val[Aiptr  ];
544                    A01=A->val[Aiptr+2];                A10=A->val[Aiptr+1];
545                    A11=A->val[Aiptr+3];                A01=A->val[Aiptr+2];
546                reg1 += A00*in1 + A01*in2;                A11=A->val[Aiptr+3];
547                reg2 += A10*in1 + A11*in2;            reg1 += A00*in1 + A01*in2;
548                }            reg2 += A10*in1 + A11*in2;
549        }        }
550        out[  2*ir] += alpha * reg1;        out[  2*irow] += alpha * reg1;
551        out[1+2*ir] += alpha * reg2;        out[1+2*irow] += alpha * reg2;
552      }      }
553      } else if (A ->col_block_size==3 && A->row_block_size ==3) {      } else if (A ->col_block_size==3 && A->row_block_size ==3) {
554          #pragma omp for private(ir,reg1,reg2,reg3,iptr,ic,Aiptr,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22) schedule(static)      int iptr0;
555      for (ir=0;ir< A->pattern->myNumOutput;ir++) {          #pragma omp for private(irow,reg1,reg2,reg3,iptr,ic,Aiptr,in1,in2,in3,A00,A10,A20,A01,A11,A21,A02,A12,A22,itmp,iptr0) schedule(static)
556        for (irow=0;irow< A->pattern->myNumOutput;irow++) {
557            reg1=0.;            reg1=0.;
558            reg2=0.;            reg2=0.;
559            reg3=0.;            reg3=0.;
560            #pragma ivdep            #pragma ivdep
561        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {        /* index[] is an increasing sequence, can find min first */
562          for (iptr0=(A->pattern->ptr[irow]);iptr0<(A->pattern->ptr[irow+1]); ++iptr0) {
563                  itmp=A->pattern->index[iptr0];
564                  if ( min_index<= itmp) { break; }
565          }
566              #pragma ivdep
567          for (iptr=iptr0;iptr<(A->pattern->ptr[irow+1]); ++iptr) {
568                itmp=A->pattern->index[iptr];                itmp=A->pattern->index[iptr];
569                if (( min_index<= itmp) && (itmp<max_index) ) {            /* index[] is an increasing sequence, can quit loop if at max */
570                ic=3*(itmp-min_index);            if (itmp >= max_index) { break; }
571                    Aiptr=iptr*9;            ic=3*(itmp-min_index);
572                    in1=in[ic];                Aiptr=iptr*9;
573                    in2=in[1+ic];                in1=in[ic];
574                    in3=in[2+ic];                in2=in[1+ic];
575                    A00=A->val[Aiptr  ];                in3=in[2+ic];
576                    A10=A->val[Aiptr+1];                A00=A->val[Aiptr  ];
577                    A20=A->val[Aiptr+2];                A10=A->val[Aiptr+1];
578                    A01=A->val[Aiptr+3];                A20=A->val[Aiptr+2];
579                    A11=A->val[Aiptr+4];                A01=A->val[Aiptr+3];
580                    A21=A->val[Aiptr+5];                A11=A->val[Aiptr+4];
581                    A02=A->val[Aiptr+6];                A21=A->val[Aiptr+5];
582                    A12=A->val[Aiptr+7];                A02=A->val[Aiptr+6];
583                    A22=A->val[Aiptr+8];                A12=A->val[Aiptr+7];
584                reg1 += A00*in1 + A01*in2 + A02*in3;                A22=A->val[Aiptr+8];
585                reg2 += A10*in1 + A11*in2 + A12*in3;            reg1 += A00*in1 + A01*in2 + A02*in3;
586                reg3 += A20*in1 + A21*in2 + A22*in3;            reg2 += A10*in1 + A11*in2 + A12*in3;
587                }            reg3 += A20*in1 + A21*in2 + A22*in3;
588        }        }
589        out[  3*ir] += alpha * reg1;        out[  3*irow] += alpha * reg1;
590        out[1+3*ir] += alpha * reg2;        out[1+3*irow] += alpha * reg2;
591        out[2+3*ir] += alpha * reg3;        out[2+3*irow] += alpha * reg3;
592      }      }
593      } else {      } else {
594          #pragma omp for private(ir,iptr,irb,icb,irow,icol,reg) schedule(static)      int iptr0;
595      for (ir=0;ir< A->pattern->myNumOutput;ir++) {          #pragma omp for private(irow,iptr,irb,icb,irtmp,ictmp,reg,itmp,iptr0) schedule(static)
596        for (iptr=A->pattern->ptr[ir];iptr<A->pattern->ptr[ir+1]; iptr++) {      for (irow=0;irow< A->pattern->myNumOutput;irow++) {
597              #pragma ivdep
598          /* index[] is an increasing sequence, can find min first */
599          for (iptr0=(A->pattern->ptr[irow]);iptr0<(A->pattern->ptr[irow+1]); ++iptr0) {
600                  itmp=A->pattern->index[iptr0];
601                  if ( min_index<= itmp) { break; }
602          }
603              #pragma ivdep
604          for (iptr=iptr0;iptr<(A->pattern->ptr[irow+1]); ++iptr) {
605                itmp=A->pattern->index[iptr];                itmp=A->pattern->index[iptr];
606                if (( min_index<= itmp) && (itmp<max_index) ) {            /* index[] is an increasing sequence, can quit loop if at max */
607               for (irb=0;irb< A->row_block_size;irb++) {            if (itmp >= max_index) { break; }
608                 irow=irb+A->row_block_size*ir;            for (irb=0;irb< A->row_block_size;irb++) {
609                     reg=0.;              irtmp=irb+A->row_block_size*irow;
610                 for (icb=0;icb< A->col_block_size;icb++) {                  reg=0.;
611               icol=icb+A->col_block_size*(itmp-min_index);              for (icb=0;icb< A->col_block_size;icb++) {
612               reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];            ictmp=icb+A->col_block_size*(itmp-min_index);
613                 }            reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[ictmp];
614                 out[irow] += alpha * reg;              }
615               }              out[irtmp] += alpha * reg;
616           }            }
617            }            }
618          }          }
619      }      }

Legend:
Removed from v.1164  
changed lines
  Added in v.1165

  ViewVC Help
Powered by ViewVC 1.1.26