symamgsol.c

a software code for computing selected eigenvalues of large sparse symmetric matrices

	 for (i=sumn; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	     buff[n+i]=buff[i]-sol[i];
#else
	     buff[n+i].r=buff[i].r-sol[i].r;
	     buff[n+i].i=buff[i].i-sol[i].i;
#endif
	 } // end for i
	 prev=next;
	 next=next->next;
	 lev++;
   } // end while




   // last level
   nB=next->nB;
   p=next->p;
   invq=next->invq;

#ifdef PRINT_INFO
   printf("%3d\n",lev);fflush(stdout);
#endif




   /* -----   last block   ----- */

   /* did we finally switch to full matrix processing? */
   if (lev>1 && next->LU.ja==NULL) {
#ifdef PRINT_INFO
      printf("dense U-factor\n");fflush(STDOUT);
      for (i=0; i<nB; i++)
	  printf("%12d",next->LU.ia[i]);
      printf("\n");fflush(stdout);
      j=0;
      for (k=0; k<nB; k++) {
	  for (i=k; i<nB; i++) {
	      printf("%12.4le",next->LU.a[j++]);
	  }
	  printf("\n");fflush(stdout);
      }
      printf("\n");fflush(stdout);
#endif

      i=1;
      COPY(&nB, buff+n+sumn,&i, sol+sumn,&i);
#ifdef PRINT_INFO
      printf("rhs\n");fflush(STDOUT);
      for (i=0; i<nB; i++)
	  printf("%12d",sumn+i+1);
      printf("\n");fflush(stdout);
      for (i=0; i<nB; i++)
	  printf("%12.4le",sol[sumn+i]);
      printf("\n");fflush(stdout);
      i=1;
#endif
      if (globalflag)
	 MYPRIVATESPTRS("L", &nB, &i, next->LU.a,next->LU.ia,
			sol+sumn, &nB, &ierr,1);     
      else
	 MYSPTRS("L", &nB, &i, next->LU.a,next->LU.ia,
		 sol+sumn, &nB, &ierr,1);     
#ifdef PRINT_INFO
      printf("sol, ierr=%d\n",ierr);fflush(STDOUT);
      for (i=0; i<nB; i++)
	  printf("%12.4le",sol[sumn+i]);
      printf("\n");fflush(stdout);
      i=1;
#endif
   }
   else {
      if (lev>1) {
	 scal=next->colscal-sumn;
	 buff+=n;
	 for (i=sumn; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	   buff[i]*=scal[i];
#else
	   val=buff[i].r;
	   buff[i].r= val*scal[i].r+buff[i].i*scal[i].i;
	   buff[i].i=-val*scal[i].i+buff[i].i*scal[i].r;
#endif
	 } // end for i
	 buff-=n;
      } /* end if */

      /* reorder the whole system */
      p-=sumn;
      for (i=sumn; i<n; i++)
	  buff[i]=buff[n+sumn+p[i]-1];
      p+=sumn;
     

      if (next->LU.ja[0]<0) {
	 flag=-1;
	 globalflag=-1;
      }
      else
	 flag=0;
      next->LU.ja[0]=IABS(next->LU.ja[0]);
      if (flag) {
	 buff+=sumn;
	 // printf("8.\n");fflush(stdout);
	 MYSYMPILUCLSOL(&nB, buff,buff+n, next->LU.a,next->LU.ja); 
	 i=0;
	 a=next->absdiag;
	 buff+=n;
	 // printf("9.\n");fflush(stdout);
	 while (i<nB) {
	       // printf("%3d\n",i);fflush(stdout);
	       if (next->LU.ja[nB+1+i]==0) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
		  // printf("%12.4le\n",buff[i]);fflush(stdout);
		  // printf("%12.4le\n",a[i]);fflush(stdout);
		  buff[i]*=a[i];
#else
		  val      =buff[i].r*a[i].r-buff[i].i*a[i].i;
		  buff[i].i=buff[i].r*a[i].i+buff[i].i*a[i].r;
		  buff[i].r=val;
#endif
		  i++;
	       }
	       else {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
		  val      =a[i]            *buff[i]+a[nB+1+i]*buff[i+1];
		  buff[i+1]=a[nB+1+i]       *buff[i]+a[i+1]   *buff[i+1];
		  buff[i]  =val;
#else
		  val        =a[i].r     *buff[i].r  -a[i].i     *buff[i].i
		             +a[nB+1+i].r*buff[i+1].r-a[nB+1+i].i*buff[i+1].i;
		  buff[i].i  =a[i].r     *buff[i].i  +a[i].i     *buff[i].r
		             +a[nB+1+i].r*buff[i+1].i+a[nB+1+i].i*buff[i+1].r;
#ifdef _COMPLEX_SYMMETRIC_		     
		  vali       =a[nB+1+i].r*buff[i].r  -a[nB+1+i].i*buff[i].i
		             +a[i+1].r   *buff[i+1].r-a[i+1].i   *buff[i+1].i;
		  buff[i+1].i=a[nB+1+i].r*buff[i].i  +a[nB+1+i].i*buff[i].r
		             +a[i+1].r   *buff[i+1].i+a[i+1].i   *buff[i+1].r;
#else
		  vali       =a[nB+1+i].r*buff[i].r  +a[nB+1+i].i*buff[i].i
		             +a[i+1].r   *buff[i+1].r-a[i+1].i   *buff[i+1].i;
		  buff[i+1].i=a[nB+1+i].r*buff[i].i  -a[nB+1+i].i*buff[i].r
		             +a[i+1].r   *buff[i+1].i+a[i+1].i   *buff[i+1].r;
#endif
		  buff[i+1].r=vali;
		  buff[i].r  =val;
#endif
		  i+=2;
	       }
	 } // end while
	 // printf("10.\n");fflush(stdout);
	 MYSYMPILUCUSOL(&nB, buff,buff, next->LU.a,next->LU.ja); 
	 // printf("11.\n");fflush(stdout);
	 buff-=n+sumn;
      }
      else
	 MYSYMPILUCSOL(&nB, buff+sumn,buff+n+sumn, next->LU.a,next->LU.ja); 
      if (flag)
	 next->LU.ja[0]=-IABS(next->LU.ja[0]);


      /* reorder the whole system back */
      buff+=n+sumn-1;
      invq-=sumn;
      for (i=sumn; i<n; i++) 
	sol[i]=buff[invq[i]];
      invq+=sumn;
      buff-=n+sumn-1;
	
      if (lev>1) {
	 scal=next->colscal-sumn;
	 for (i=sumn; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	     sol[i]*=scal[i];
#else
	     val=sol[i].r;
	     sol[i].r=val*scal[i].r-sol[i].i*scal[i].i;
	     sol[i].i=val*scal[i].i+sol[i].i*scal[i].r;
#endif
	 } // end for i
      } /* end if */
      
   } // end if-else (next->LU.ja==NULL)

      

   /* block backward substitution */
   lev=nlev-1;
   while (lev>0) {
	 nB=prev->nB;
	 invq=prev->invq;

#ifdef PRINT_INFO
	 printf("%3d\n",lev);fflush(stdout);
#endif
	 
	 /* copy r.h.s */
	 for (i=sumn; i<n; i++)
	     buff[i]=sol[i];

	 /* second block */
	 MYCSRMATVEC(prev->F,sol+sumn,buff+sumn-nB);

	 /* first block */
	 sumn-=nB;


	 if (prev->LU.ja[0]<0) {
	    flag=-1;
	    globalflag=-1;
	 }
	 else
	    flag=0;
	 prev->LU.ja[0]=IABS(prev->LU.ja[0]);
	 if (flag) {
	    sol +=sumn;
	    buff+=sumn;
	    // printf("12.\n");fflush(stdout);
	    MYSYMPILUCLSOL(&nB, buff,buff, prev->LU.a,prev->LU.ja);
	    // printf("13.\n");fflush(stdout);
	    i=0;
	    a=prev->LU.a;
	    while (i<nB) {
	          if (prev->LU.ja[nB+1+i]==0) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
		     sol[i]=buff[i]*a[i];
#else
		     sol[i].r=buff[i].r*a[i].r-buff[i].i*a[i].i;
		     sol[i].i=buff[i].r*a[i].i+buff[i].i*a[i].r;
#endif
		     i++;
		  }
		  else {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
		     sol[i]  =a[i]            *buff[i]+a[nB+1+i]*buff[i+1];
		     sol[i+1]=a[nB+1+i]       *buff[i]+a[i+1]   *buff[i+1];
#else
		     sol[i].r  =a[i].r     *buff[i].r  -a[i].i     *buff[i].i
		               +a[nB+1+i].r*buff[i+1].r-a[nB+1+i].i*buff[i+1].i;
		     sol[i].i  =a[i].r     *buff[i].i  +a[i].i     *buff[i].r
		               +a[nB+1+i].r*buff[i+1].i+a[nB+1+i].i*buff[i+1].r;
#ifdef _COMPLEX_SYMMETRIC_		     
		     sol[i+1].r=a[nB+1+i].r*buff[i].r  -a[nB+1+i].i*buff[i].i
		               +a[i+1].r   *buff[i+1].r-a[i+1].i   *buff[i+1].i;
		     sol[i+1].i=a[nB+1+i].r*buff[i].i  +a[nB+1+i].i*buff[i].r
		               +a[i+1].r   *buff[i+1].i+a[i+1].i   *buff[i+1].r;
#else
		     sol[i+1].r=a[nB+1+i].r*buff[i].r  +a[nB+1+i].i*buff[i].i
		               +a[i+1].r   *buff[i+1].r-a[i+1].i   *buff[i+1].i;
		     sol[i+1].i=a[nB+1+i].r*buff[i].i  -a[nB+1+i].i*buff[i].r
		               +a[i+1].r   *buff[i+1].i+a[i+1].i   *buff[i+1].r;
#endif
#endif
		     i+=2;
		  }
	    } // end while
	    // printf("14.\n");fflush(stdout);
	    sol -=sumn;
	    buff-=sumn;
	 }
	 else
	    MYSYMPILUCSOL(&nB, buff+sumn,sol+sumn, prev->LU.a,prev->LU.ja);

	 /* update solution */ 
	 buff+=sumn;
	 sol+=sumn;
	 for (i=0; i<nB; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	     buff[i]=buff[n+i]-sol[i];
#else
	     buff[i].r=buff[n+i].r-sol[i].r;
	     buff[i].i=buff[n+i].i-sol[i].i;
#endif
	 } // end for i
	 if (flag) {
	    // printf("15.\n");fflush(stdout);
	    MYSYMPILUCUSOL(&nB, buff,buff, prev->LU.a,prev->LU.ja);
	    // printf("16.\n");fflush(stdout);
	    prev->LU.ja[0]=-IABS(prev->LU.ja[0]);
	 }
	 sol-=sumn;
	 buff-=sumn;

	 /* reorder the whole system */
	 buff+=sumn-1;
	 invq-=sumn;
	 for (i=sumn; i<n; i++)
	     sol[i]=buff[invq[i]];
	 invq+=sumn;
	 buff-=sumn-1;

	 if (lev>1) {
	    scal=prev->colscal-sumn;
	    for (i=sumn; i<n; i++) {
#if defined _DOUBLE_REAL_ || defined _SINGLE_REAL_
	        sol[i]*=scal[i];
#else
		val=sol[i].r;
		sol[i].r=val*scal[i].r-sol[i].i*scal[i].i;
		sol[i].i=val*scal[i].i+sol[i].i*scal[i].r;
#endif
	    } // end for i
	 } /* end if */

	 
	 next=prev;
	 prev=prev->prev;
	 lev--;
   } // end while


} /* end SYMAMGsol */